text:
function img_pad=ImgPad(img_in, pad_size, op_type, varargin)
%Perform matrix border padding and the corresponding inverse transform
% Input: mat_in - the input matrix
%             pad_size - the size of padding
%             op_type - the operation type (0:padding, 1: recovering)
%             varargin: the padding method when op_type=0
%Output: mat_out - result of the operation
% 2016-10-19 jlfeng
[nr,nc]=size(img_in);
pad_size_half=ceil(pad_size/2);
if (0==op_type)
    switch(varargin{1})
        case 0
            pad_type='symmetric';
        case 1
            pad_type='replicate';
        case 2
            pad_type='replicate';
    end
    img_pad=padarray(img_in,[pad_size_half,pad_size_half],pad_type,'both');
elseif (1==op_type)
    img_pad=img_in(pad_size_half+1:nr-pad_size_half,pad_size_half+1:nc-pad_size_half);
else
    error('Unknown operation type.');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function rs = diff(s, nth)

%tstoolbox/@signal/diff
%   Syntax:
%     * diff(s, nth)
%
%   Compute the nth numerical derivative along dimension 1. s has be to
%   sampled equidistantly.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,2);

if nargin < 2
	nth = 1;
end

a = getaxis(s,1);
switch  resolution(a)
	case 'linear'
		c = diff(s.core, nth, delta(a)); 		% call real working routine for parent core object
		rs = signal(c, s);				% special constructor calling syntax for working routines
		rs = setyunit(rs, yunit(s)/(unit(a)^nth));
		a = setfirst(a, first(a) + nth*delta(a)/2);
		rs = setaxis(rs,  1, a);
		rs = addhistory(rs,  ['Computed ' num2str(nth) '# numerical derivative along dimension 1'] );
		rs = addcommandlines(rs, 's = diff(s', nth);

        otherwise
	   x=spacing(s)';
	   x2=x(1:end-1);
	   x=x(2:end);

	   y=data(s);
	   y2=y(1:end-1,:);
	   y=y(2:end,:);
	   y=(y-y2);
	   for i=1:length(y(1,:))
	     y(:,i)=y(:,i)./(x-x2);
	   end
	   x=(x-x2)/2+x2;

	   rs=signal(core(y),s);
	   rs=setplothint(rs,'multigraph');
	   rs = setyname(rs, 'd ld N(r)/d lg r');
	   rs = setyunit(rs, yunit(s)/(unit(a)^nth));
	   rs = setaxis(rs,  1, achse(x));
	   rs = addhistory(rs,  ['Computed ' num2str(nth) '# numerical derivative along dimension 1'] );
	   rs = addcommandlines(rs, 's = diff(s', nth);

%	   error('Data values are not sampled equidistantly');

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function b = uniquestripped(a)
%UNIQUESTRIPPED  Internal function (version without checkings etc.)

b = sort(a(:)');
i = diff([b NaN])~=0;
if ~all(i)
    b = b(i);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function x = triu( x, k )

%   Disciplined convex/geometric programming information for TRIU:
%       TRIU imposes no convexity restrictions on its arguments.

%
% Check inputs
%

if nargin < 2, k = 0; end
s = x.size_;
if length( s ) > 2,
    error( 'The first argument must be 2-D.' );
elseif ~isnumeric( k ) || length( k ) ~= 1,
    error( 'The second argument must be an integer scalar.' );
end

%
% Zero out the elements outside of the desired triangle
%

b = x.basis_;
b( :, ~triu(ones(s),k) ) = 0;
x = cvx( s, b );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function net=train_denoising_mini(source_1, source_2, context_win, hidden_units, num_layers, isdropout, ...
    isRNN, iscleanonly, circular_step , isinputL1, MFCCorlogMelorSpectrum, ...
    framerate, pos_neg_r, outputnonlinear, opt, act, train_mode, const,  ...
    const2, isGPU, max_iter, batchsize, lbfgs_iter, grad_clip, wdecay)
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
% Demo Denoising training ---------------------------------------------
% context_win - context window size
% hidden_units - hidden units
% num_layers - layer number
% isdropout - 1: use dropout, 0: no dropout
% isRNN - RNN temporal connection
% iscleanonly - One output source or two
% circular_step - Circular shift step
% isinputL1 - normalize input as L1 norm = 1
% MFCCorlogMelorSpectrum - 0: MFCC, 1: logmel, 2: spectra
% framerate - feature frame rate
% pos_neg_r - discriminative training gamma parameter
% outputnonlinear - Last layer - linear or nonlinear
% softabs - soft mask obj
% act - 0: logistic, 1: tanh, 2: RELU
% const - constant for avoiding numerical problems
% const2- constant for avoiding numerical problems
% isGPU - 0: not using GPU, 1: using GPU
% train_mode - 0
% max_iter - max LBFGS iterations

rand('state',0)
randn('state',0)

%% setup paths for code. assumed this script runs in its own directory
% CHANGE baseDir to the top of code directory
baseDir= ['..',filesep, '..',  filesep];
codeDir = [baseDir,'codes', filesep];
minFuncDir = [baseDir, 'tools', filesep, 'minFunc_2012', filesep];

saveDir = [codeDir,filesep,'denoising',filesep,'discrim_joint_offset_results'];

%% add paths
addpath([baseDir, filesep,'tools', filesep,'labrosa']);
% addpath([baseDir, filesep,'tools', filesep,'bss_eval']);
addpath([baseDir, filesep,'tools', filesep,'bss_eval_2.1']);
% addpath([baseDir, filesep,'tools', filesep,'bss_eval_3']);
addpath(baseDir);
addpath(genpath(minFuncDir));
addpath(codeDir);
addpath([codeDir,'denoising']);

CFGPath=[baseDir,'tools',filesep,'htk_features', filesep];
addpath(CFGPath);

%% setup network architecture
setup_architecture;

%% initialize weights
[stack_i, W_t_i] = initialize_weights(eI);
[theta] = rnn_stack2params(stack_i, eI, W_t_i);

%% Directory of features
eI.featInBase =baseDir;

%% load data
eI.useCache = 0;

%% setup minFunc
options.Diagnostics = 'on';
options.Display = 'iter';
% options.MaxIter = 400;

options.MaxIter = lbfgs_iter;
options.MaxFunEvals = 2500;
options.Corr = 50;
options.DerivativeCheck = 'off';
% options.DerivativeCheck = 'on';
options.outputFcn = @save_callback_denoising_general;

% eI.DataPath=[codeDir,'mir1k', filesep, 'Wavfile',filesep];
eI.iterStart=-batchsize;

%% feed data
for ii= 1: max_iter
    fprintf('Global iterations: %d\n', ii);
    eI.iterStart=eI.iterStart+batchsize;

    batch_cell_idx = randperm(numel(source_1), min(batchsize, numel(source_1)));
    source_1_mini = cell(1, numel(batch_cell_idx)); source_2_mini = cell(1, numel(batch_cell_idx));
    for jj=1:numel(batch_cell_idx)
        source_1_mini{jj} = source_1{batch_cell_idx(jj)};
        source_2_mini{jj} = source_2{batch_cell_idx(jj)};
    end
    % train_files= dir( [eI.DataPath, 'train',filesep,'*wav']);
    % 0 -- chunk, 2--no chunk
    [data_cell, targets_cell, mixture_spectrum] = ...
                formulate_data(source_1_mini, source_2_mini, eI, eI.train_mode);

    %% BSS EVAL setting
    eI.writewav=0;

    %% run optimizer
    if isfield(eI,'cleanonly') && eI.cleanonly==1,
        % % for non joint training
        eI.isdiscrim=1;
        [theta,val]=minFunc(@drdae_discrim_obj, theta, options, eI, data_cell, ...
                            targets_cell, false, false);
    else
       if isGPU==1 && strcmpi(eI.activationFn,'RELU') && opt==1
          [theta,val]=minFunc(@drdae_discrim_joint_kl_obj_gpu_relu, theta, options, eI, ...
           data_cell, targets_cell, mixture_spectrum, false, false);
       elseif isGPU==1
          [theta,val]=minFunc(@drdae_discrim_joint_kl_obj_gpu, theta, options, eI, ...
            data_cell, targets_cell, mixture_spectrum, false, false);
       else
          [theta,val]=minFunc(@drdae_discrim_joint_kl_obj, theta, options, eI, ...
            data_cell, targets_cell, mixture_spectrum, false, false);
       end
    end

end

net.theta = theta;
net.eI = eI;

return;

%% unit test - small example
% context window size
context_win = 1;
% hidden units
hidden_units = 16;
num_layers = 1;
isdropout = 0;
% RNN temporal connection
isRNN = 2;
% One output source or two
iscleanonly = 0;
% Circular shift step
circular_step = 1000000;
% normalize input as L1 norm = 1
isinputL1 = 0;
% 0: MFCC, 1: logmel, 2: spectra
MFCCorlogMelorSpectrum = 2;
% feature frame rate
framerate = 64;

% discriminative training gamma parameter
pos_neg_r = 0.05;
% Last layer - linear or nonlinear
outputnonlinear = 0;
% soft mask obj
softabs = 1;
% 0: logistic, 1: tanh, 2: RELU
act = 2;
% constant for avoiding numerical problems
const = 1e-10;
% constant for avoiding numerical problems
const2 = 0.001;
% 0: not using GPU, 1: using GPU
isGPU = 0;

train_mode = 0;
% 0:'softlinear',1:'softabs', 2:'softquad', 3:'softabs_const',
% 4:'softabs_kl_const'
opt = 1;

[train1, fs, nbits]=wavread('female_train.wav');
[train2, fs, nbits]=wavread('male_train.wav');
maxLength=max([length(train1), length(train2)]);
train1(end+1:maxLength)=eps;
train2(end+1:maxLength)=eps;

max_iter = 30;

train_denoising({train1'}, {train2'}, context_win, hidden_units, num_layers, isdropout, ...
    isRNN, iscleanonly, circular_step , isinputL1, MFCCorlogMelorSpectrum, ...
    framerate, pos_neg_r, outputnonlinear, opt, act, train_mode, const,  ...
    const2, isGPU, max_iter)

%% unit test 2 - best setting:
% context window size
context_win = 3;
% hidden units
hidden_units = 1000;
num_layers = 3;
isdropout = 0;
% RNN temporal connection
isRNN = 2;
% One output source or two
iscleanonly = 0;
% Circular shift step
circular_step = 10000;
% normalize input as L1 norm = 1
isinputL1 = 0;
% 0: MFCC, 1: logmel, 2: spectra
MFCCorlogMelorSpectrum = 2;
% feature frame rate
framerate = 64;
% discriminative training gamma parameter
pos_neg_r = 0.05;
% Last layer - linear or nonlinear
outputnonlinear = 0;
% soft mask obj
softabs = 1;
% 0: logistic, 1: tanh, 2: RELU
act = 2;
% constant for avoiding numerical problems
const = 1e-10;
% constant for avoiding numerical problems
const2 = 0.001;
% 0: not using GPU, 1: using GPU
isGPU = 0;

train_mode = 0;
% 0:'softlinear',1:'softabs', 2:'softquad', 3:'softabs_const',
% 4:'softabs_kl_const'
opt = 1;

[train1, fs, nbits]=wavread('female_train.wav');
[train2, fs, nbits]=wavread('male_train.wav');
maxLength=max([length(train1), length(train2)]);
train1(end+1:maxLength)=eps;
train2(end+1:maxLength)=eps;

max_iter = 30;

train_denoising({train1'}, {train2'}, context_win, hidden_units, num_layers, isdropout, ...
    isRNN, iscleanonly, circular_step , isinputL1, MFCCorlogMelorSpectrum, ...
    framerate, pos_neg_r, outputnonlinear, opt, act, train_mode, const,  ...
    const2, isGPU, max_iter)

author: posenhuang

repo: deeplearningsourceseparation


text:
function [out] = phenology_2(p1,p2,p3,t,tmax,dt)
%phenology_2

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Phenology-based maximum interception capacity (returns store size [mm])
% Constraints:  Implicit assumption: 0 <= p2 <= 1
% @(Inputs):    p1   - mean interception capacity [mm]
%               p2   - seasonal change as fraction of the mean [-]
%               p3   - time of maximum store size [d]
%               t    - current time step [-]
%               tmax - seasonal length [d]
%               dt   - time step size [d]

out = p1*(1+p2*sin(2*pi*(t*dt-p3)/tmax));

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function m=mapzo(x)
xmin=min(x);xmax=max(x);
for i=1:length(x)
    m(i)=(x(i)-xmin)/(xmax-xmin);
end

author: bastamon


text:
function f = p40_f ( n, x )

%*****************************************************************************80
%
%% P40_F evaluates the objective function for problem 40.
%
%  Discussion:
%
%    There is a typo in the reference.  I'm just guessing at the correction.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  f = x(1)^2 + 2.0 * x(2)^2 ...
    - 0.3 * cos ( 3.0 * pi * x(1) ) ...
    + cos ( 4.0 * pi * x(2) ) + 0.3;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = process_arima( varargin )
% PROCESS_ARIMA: Auto-regressive Moving Average filter
%
% USAGE:   sProcess = process_arima('GetDescription')
%            sInput = process_arima('Run', sProcess, sInput, method=[])
%                 F = process_arima('Compute', F, Fbase, Order)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'ARIMA filter';
    sProcess.FileTag     = 'arima';
    sProcess.Category    = 'Filter';
    sProcess.SubGroup    = 'Pre-process';
    sProcess.Index       = 68;
    sProcess.Description = '';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'raw', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'raw', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.processDim  = 1;   % Process channel by channel

    % Definition of the options
    % === Baseline time window
    sProcess.options.baseline.Comment = 'Baseline:';
    sProcess.options.baseline.Type    = 'baseline';
    sProcess.options.baseline.Value   = [];
    % === Filter order
    sProcess.options.order.Comment = 'Order of the filter:';
    sProcess.options.order.Type    = 'value';
    sProcess.options.order.Value   = {5,'',0};
    % === Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    sProcess.options.sensortypes.InputTypes = {'data', 'raw'};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Order = sProcess.options.order.Value{1};
    if isfield(sProcess.options, 'baseline') && isfield(sProcess.options.baseline, 'Value') && iscell(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value{1})
        Time = sProcess.options.baseline.Value{1};
    else
        Time = [];
    end
    if isempty(Time)
        Comment = 'ARIMA(%d,1,0) - Baseline: [All file]';
    elseif any(abs(Time) > 2)
        Comment = sprintf('ARIMA(%d,1,0) - Baseline: [%1.3fs,%1.3fs]', Order, Time(1), Time(2));
    else
        Comment = sprintf('ARIMA(%d,1,0) - Baseline: [%dms,%dms]', Order, round(Time(1)*1000), round(Time(2)*1000));
    end
end

%% ===== RUN =====
function sInput = Run(sProcess, sInput) %#ok<DEFNU>
    % Get options
    Order = sProcess.options.order.Value{1};
    if isfield(sProcess.options, 'baseline') && isfield(sProcess.options.baseline, 'Value') && iscell(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value) && ~isempty(sProcess.options.baseline.Value{1})
        BaselineBounds = sProcess.options.baseline.Value{1};
    else
        BaselineBounds = [];
    end
    % Get baseline indices
    if ~isempty(BaselineBounds)
        iBaseline = panel_time('GetTimeIndices', sInput.TimeVector, BaselineBounds);
        if isempty(iBaseline)
            bst_report('Error', sProcess, [], 'Invalid baseline definition.');
            sInput = [];
            return;
        end
    % Get all file
    else
        iBaseline = 1:size(sInput.A,2);
    end
    % Filter data
    sInput.A = Compute(sInput.A, sInput.A(:,iBaseline), Order);
    % Error handling
    if isempty(sInput.A)
        bst_report('Error', sProcess, [], 'Error while filtering the signal.');
        sInput = [];
        return;
    end
    % Do not keep the Std field in the output
    if isfield(sInput, 'Std') && ~isempty(sInput.Std)
        sInput.Std = [];
    end
end

%% ===== EXTERNAL CALL =====
% USAGE: process_arima('Compute', F, Fbase=[], Order=5)
function F = Compute(F, Fbase, Order)
    % Default order: 5
    if (nargin < 3) || isempty(Order)
        Order = 5;
    end
    % If there is no baseline, use the whole time segment
    if (nargin < 2) || isempty(Fbase)
        Fbase = F;
    end
    Nsig = size(F,1);

    % If order is 0, just diff the signal
    if (Order == 0)
        % Detrend and diff data to filter
        F = diff(detrend(F'))';
        % Add one sample at the beginning to account for the diff
        F = [F(:,2), F];
    else
        % Detrend and diff baseline
        Fbase = diff(detrend(Fbase'));
        % Compute AR model for each signal
        if bst_get('UseSigProcToolbox')
            arm = lpc(Fbase, Order);
        else
            arm = zeros(Nsig, Order+1);
            for i = 1:Nsig
                arm(i,:) = oc_lpc(Fbase(:,i), Order);
            end
        end
        % Remove the ones for each it cannot be estimated
        arm(any(isnan(arm),2),:) = [];
        % If there is nothing left, error
        if isempty(arm)
            F = [];
        end
        % Average all the models
        arm = mean(arm,1);

        % Detrend and diff data to filter
        F = diff(detrend(F'));
        % Add a few samples of mirrored signal at the beginning to minimize the edge effects
        F = [F(length(arm)+1:-1:1,:); F];
        % Apply filter to the data
        F = filter(arm, 1, F)';
        % Remove the mirrored samples (keep one, to account for the diff)
        F(:,1:length(arm)) = [];
    end
end

author: brainstorm-tools


text:
function alpha = get_Clambda_zero_points(v)
alpha = zeros(v + 1, 1);
alpha(1) = 1;
alpha(v + 1) = realmax;
%above two values are obtained by simple calculations and avoid numerical
%problems
%Newton descend method
epsilon = 1e-6;%tolerance error
for i = 2 : v
    beta = (i - 1)/v;
    x0 = 0;
    x1 = 1.5;%initial point, the zero point of d^2C(lambda)/x^2
    while(abs(x0 - x1) > epsilon)
        x0 = x1;
        x1 = x0 - (Clambda(x0) - beta)/dClambdadx(x0);
    end
    alpha(i) = x1;
end
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

h = sigmoid(X * theta);

J = (1 / m) * sum( (-y .* log(h)) - ((1 - y) .* log( 1 - h )) );

grad = (1 / m) * X' * (h - y);

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

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
%compare max, min a mean manipulability for vanilla
load('/Users/arturogilaparicio/Desktop/arte/robots/RETHINK/sgo_v0.5/experiment1/mat/experiment1A.mat')

[correct_manips]=eval_experiment(Gout, random_manips);

[y,i]=max(sum(correct_manips'));
max_manip = correct_manips(i,:);

[y,i]=min(sum(correct_manips'));
min_manip = correct_manips(i,:);

mean_manip = mean(correct_manips);

figure, plot(max_manip), hold
plot(min_manip)
plot(mean_manip)

repo: ARTE


text:
function g = rbfhKernGradient(kern, x1, x2, covGrad)

% RBFHKERNGRADIENT Gradient of RBFH kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% radial basis function heat kernel's parameters. As well as the kernel
% structure and the input positions, the user provides a matrix PARTIAL
% which gives the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations for which the gradients are being
% computed.
% ARG covGrad : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG covGrad : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO: rbfhKernParamInit, kernGradient
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 4
    covGrad = x2;
    x2 = x1;
end
if size(x1, 2) ~= 2 || size(x2, 2) ~= 2
    error('Input can only have two columns');
end

% Split the domain into time domain and spatial domain and account for
% missing values. If there are no missing values the computation of the
% kernel is a pointwise prodruct, otherwise it is a kronecker product.
t1 = x1(x1(:,1)~=Inf,1);
t2 = x2(x2(:,1)~=Inf,1);
s1 = x1(x1(:,2)~=Inf,2);
s2 = x2(x2(:,2)~=Inf,2);
if (length(t1) == length(s1)) && (length(t2) == length(s2))
    ut1 = unique(t1);
    ut2 = unique(t2);
    us1 = unique(s1);
    us2 = unique(s2);
    if (length(ut1)*length(us1) == length(t1)) && ...
            (length(ut2)*length(us2) == length(t2))
        t1 = ut1; s1 = us1; t2 = ut2; s2 = us2;
        isPointwise = false;
    else
        isPointwise = true;
    end
else
    isPointwise = false;
end

kern.rbf.inverseWidth = kern.inverseWidthTime;
Kt = rbfKernCompute(kern.rbf, t1, t2);
kern.rbf.inverseWidth = kern.inverseWidthSpace;
Ks = rbfKernCompute(kern.rbf, s1, s2);

if isPointwise
    covGradt = covGrad.*Ks;
    covGrads = covGrad.*Kt;
else
    covGradt = zeros(length(t1), length(t2));
    covGrads = zeros(length(s1), length(s2));
    endOne = 0;
    startOne = 1;
    for k=1:length(t1)
        endOne = endOne + length(s1);
        startTwo = 1;
        endTwo = 0;
        for l=1:length(t2)
            endTwo = endTwo + length(s2);
            covGradt(k,l) = sum(sum(covGrad(startOne:endOne, startTwo:endTwo).*Ks));
            startTwo = endTwo + 1;
        end
        startOne = endOne + 1;
    end
    for k=1:length(s1)
        indRows = (k:length(s1):(length(t1)*length(s1)))';
        for l=1:length(s2)
            indCols = l:length(s2):(length(t2)*length(s2));
            covGrads(k,l) = sum(sum(covGrad(indRows, indCols).*Kt));
        end
    end
end

kern.rbf.inverseWidth = kern.inverseWidthTime;
gt = rbfKernGradient(kern.rbf, t1, t2, covGradt);
kern.rbf.inverseWidth = kern.inverseWidthSpace;
gs = rbfKernGradient(kern.rbf, s1, s2, covGrads);

g = [gt(1) gs(1)];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [K, Kbase, n2]  = gaussianKernCompute(kern, x, x2)

% GAUSSIANKERNCOMPUTE Compute the Gaussian kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the Gaussian kernel given
% inputs associated with rows and columns.
% RETURN K : the kernel matrix computed at the given points.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG X : the input matrix associated with the rows of the kernel.
% ARG X2 : the input matrix associated with the columns of the kernel.
%
% FORMAT
% DESC computes the kernel matrix for the Gaussian kernel given a design matrix of inputs.
% RETURN K : the kernel matrix computed at the given points.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
%
% SEEALSO : gaussianKernParamInit, kernCompute, kernCreate, gaussianKernDiagCompute
%
% COPYRIGHT : Mauricio Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS: Mauricio Alvarez, 2009

% KERN

if kern.isArd
    sqrtP = sqrt(kern.precisionU);
    sqrtPx = x*sparseDiag(sqrtP);
    if nargin < 3
        n2 = dist2(sqrtPx, sqrtPx);
    else
        sqrtPx2 = x2*sparseDiag(sqrtP);
        n2 = dist2(sqrtPx, sqrtPx2);
    end
    Kbase = exp(-0.5*n2);
else
    if nargin < 3
        n2 = dist2(x, x);
    else
        n2 = dist2(x, x2);
    end
    Kbase = exp(-0.5*kern.precisionU*n2);
end
K = kern.sigma2Latent*Kbase;


text:
function [nodeInterface, elemInterface, nodeBg, elemBg] = get2dpolymesh(node, elem, phi)
%% GET2DPOLYMESH gets the interface fitted polygonal mesh for a 2d quadtree mesh

%% Prepare
if ~iscell(elem); elem = num2cell(elem,2); end

elemBg = elem;
nodeBg = node;
T = auxstructurepoly(elem);
edge = T.edge;
edge2elem = double(T.edge2elem);

N = size(node,1);
NT = size(elem,1);

%% compute the level set value at each vertex
%
phiValue = phi(node);
vSign = msign(phiValue);

%% Step 1: Find points ON interface
% Find the intersection points between edges and the interface
isCutEdge = (vSign(edge(:,1)).*vSign(edge(:,2))<0);
A = node(edge(isCutEdge,1),:);
B = node(edge(isCutEdge,2),:);
nodeCut = findintersectbisect(phi,A,B);
Ncut = size(nodeCut, 1);
ixCutNode = N+(1:Ncut)';
ixIntersectNode = find(abs(vSign) <= eps);
Nintersect = size(ixIntersectNode,1);
vSign(N+1:N+Ncut) = 0;
vSign(N+Ncut+1:N+Ncut+Nintersect) = 0;

%% Step 2: find interface elem and nodes
isInterfaceElem = false(NT,1);
interfaceElem = edge2elem(isCutEdge,[1,2]);
isInterfaceElem(interfaceElem) = true;
% new vertices needed to be merged into interface element
elemUpdate = accumarray(interfaceElem(:),[ixCutNode;ixCutNode], [NT,1],@(x){x'});
elem = cellfun(@horzcat, elem, elemUpdate, 'UniformOutput', false);
elemVertNum = cellfun('length',elem);

vSign2elem = mat2cell(vSign([elem{:}]),elemVertNum);
vSign2elem = cellfun(@transpose, vSign2elem, 'UniformOutput', 0);
isInterfaceElem(cellfun(@(x) sum(x==0)==2, vSign2elem)) = true; % 2 vertices on interface
ixIFElem = find(isInterfaceElem);
% elemOnInterface = elem(isInterfaceElem,:);

nodeInterface = [node; nodeCut];

%% cut polygonal element into sub elements for computation
numInterfaceElem = length(ixIFElem);
elem(NT+1:NT+numInterfaceElem) = cell(numInterfaceElem,1);

for i = 1:numInterfaceElem
   ix = ixIFElem(i);
   ixLocElem1= [find(vSign2elem{ix}==0), find(vSign2elem{ix}==1)];
   ixLocElem2 = [find(vSign2elem{ix}==0), find(vSign2elem{ix}==-1)];
%    elem{ix} = {elem{ix}(ixLocElem1); elem{ix}(ixLocElem2)};
    elem{NT+i} = elem{ix}(ixLocElem2);
    elem{ix} = elem{ix}(ixLocElem1);
end

% ixIFElemNew = [ixIFElem; ((NT+1):(NT+numInterfaceElem))'];
elemInterface = fixorientationpoly(nodeInterface,elem);

end

repo: ifem


text:
function displaySegmentGraph(im, segimage, adjlist)

stats = regionprops(segimage, 'Centroid');
centroids = cat(1, stats.Centroid);

segimage = double(segimage) / double(max(segimage(:)));

[gx, gy] = gradient(segimage);

edges = find((gx ~= 0) | (gy ~= 0));

[imh, imw, imb] = size(im);

for b = 1:imb
    im(edges+(b-1)*imh*imw) = 1*(b==1); % make edges red
end

hold off, imagesc(im), axis image, hold on

% convert adjacency matrix to adjacency list if necessary
if size(adjlist, 1) == size(adjlist, 2)
    adjmat = adjlist;
    [s1, s2] = find(adjmat);
    adjlist = [s1 s2];
    ind = find(adjlist(:, 1) < adjlist(:, 2));
    adjlist = adjlist(ind, :);
end

if size(adjlist, 2) == 2
    for k = 1:size(adjlist, 1)
        plot(centroids(adjlist(k, :), 1), centroids(adjlist(k, :), 2), 'g', 'LineWidth', 1);
    end
end

author: zouchuhang

repo: LayoutNet


text:
function [p] = mesh_scalp_interp(p)

% mesh_scalp_interp - Interpolates scalar over scalp mesh from electrodes
%
% Usage: [p] = mesh_scalp_interp(p)
%
% p is the eeg_toolbox struct (see eeg_toolbox_defaults)
%
% This function returns the laplacian matrix of the scalp
% mesh, the interpolation matrix based on the min norm of
% the laplacian, and the voltage interpolated on the scalp.
% It also calls an electrode "coregistration" function and
% the results of that operation.
%
% The return matrices are all in p.mesh.data, including the
% interpolated timeseries values in p.mesh.data.timeseries{N},
% where N is the index of the scalp mesh in this case (see
% mesh_check to get N).
%
% This function is called by eeg_contours_engine.
%
% See also, mesh_fit_elec, mesh_laplacian, mesh_laplacian_interp
% for more information.
%

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

% Licence:  Gnu GPL, no express or implied warranties
% History:  09/02, Darren.Weber_at_radiology.ucsf.edu
%                  extracted out of MESH_SCALP_INTERP, called by that
%                  function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('p','var'),
    fprintf('Setting default parameters.\n');
   [p] = eeg_toolbox_defaults('create');
end

% Calculate voltage at scalp mesh vertices
if isempty(p.mesh.data),
    if ~isempty(p.mesh.file),
       [p] = mesh_open(p);
    else
        msg = sprintf('MESH_SCALP_INTERP: No meshes, load meshes first.');
        warning(msg);
    end
end

% get the scalp mesh
[p.mesh.current,meshExists] = mesh_check(p,'scalp');

if isempty(p.mesh.current),
    msg = sprintf('MESH_SCALP_INTERP: No ''scalp'' mesh, check p.mesh.data.meshtype.');
    error(msg);
end

% Find nearest vertices of scalp mesh to electrodes
% This will also reorder the scalp mesh vertices/faces
% so that the first 1:Nelec vertices are the vertices
% in the scalp that lie closest to the electrodes.
% Serious problems arise if there are any duplicate
% scalp vertices for any electrode, but this routine
% will continue regardless, at present.
p = mesh_fit_elec(p);

% Ensure scalp is the currently selected mesh
elecN = mesh_check(p,'elec');
scalpN = mesh_check(p,'scalp');
[p.mesh.current,meshExists] = mesh_check(p,'scalp');

% Get vertices and faces of scalp mesh, after they
% have been rearraned by mesh_fit_elec
scalpvert = p.mesh.data.vertices{scalpN};
scalpface = p.mesh.data.faces{scalpN};

% Get the scalp vertices and faces
% that correspond to the electrodes
elecvert  = p.mesh.data.vertices{elecN};
elecface  = p.mesh.data.faces{elecN};
% Get the index of each electrode vertex
% into the scalp mesh, should be first
% 1:Nelec indices of scalp vertices
elecindex = p.mesh.data.elecindex{scalpN};

fprintf('MESH_SCALP_INTERP: Calculating Scalp Interpolation\n');

% Could identify the elements of the scalp below
% the electrodes and discard them before calculating
% the interpolation matrices, or use them in the
% interpolation, but set their potential to zero
% later (see comment below).

% Calculate the Laplacian matrix of the scalp mesh (once off calculation)
if isfield(p.mesh.data,'lapmat'),
    if isempty(p.mesh.data.lapmat{scalpN}),
        p.mesh.data.lapmat{scalpN} = mesh_laplacian(scalpvert,scalpface);
    end
else
    p.mesh.data.lapmat{scalpN} = mesh_laplacian(scalpvert,scalpface);
end

% Calculate the electrode to scalp mesh interpolation matrix
% using a Laplacian minimum norm constraint (once off calculation)
if isfield(p.mesh.data,'lapint'),
    if isempty(p.mesh.data.lapint{scalpN}),
        [p.mesh.data.lapint{scalpN}, p.mesh.data.keepindex{scalpN}, p.mesh.data.repindex{scalpN}] = mesh_laplacian_interp(p.mesh.data.lapmat{scalpN}, elecindex);
    end
else
    [p.mesh.data.lapint{scalpN}, p.mesh.data.keepindex{scalpN}, p.mesh.data.repindex{scalpN}] = mesh_laplacian_interp(p.mesh.data.lapmat{scalpN}, elecindex);
end

fprintf('MESH_SCALP_INTERP: Multiplying Voltage by Interpolation matrix...\n');
tic;
% Interpolate voltage to all scalp mesh vertices (once off calc)
if isempty(p.mesh.data.repindex{scalpN}),
    p.mesh.data.Cdata{scalpN} = p.mesh.data.lapint{scalpN} * p.volt.data';
    %p.mesh.data.Cdata{scalpN} = p.volt.data * p.mesh.data.lapint{scalpN}';
else
    p.mesh.data.Cdata{scalpN} = p.mesh.data.lapint{scalpN} * p.volt.data(:,p.mesh.data.keepindex{scalpN})';
    %p.mesh.data.Cdata{scalpN} = p.volt.data(:,p.mesh.data.keepindex{scalpN}) * p.mesh.data.lapint{scalpN}';
end
t = toc;
fprintf('done (%6.2f sec).\n',t);

% Now a quick test of the interpolation: the values at the first
% Nelec vertices of the scalp mesh must be equal to the values
% of the original voltage data
TMP = p.mesh.data.Cdata{scalpN};
Nelec = size(p.volt.data,2);
if ~isequal(p.volt.data',TMP(1:Nelec,:)),
    msg = sprintf('MESH_SCALP_INTERP: Fatal Error in Mesh Voltage Interpolation, Reload all Data?\n');
    error(msg);
end

% At present, all scalp vertices are interpolated.  It might
% be useful at this point to determine all vertices that are
% "well below" the electrode array, say > 5cm from a nearest
% electrode and then set their voltage to zero.

return

author: PatternRecognition

repo: OpenBMI


text:
%--------------------------------------------------------------------------------------------------------
% The system is created based on the principles described in the following paper
% Jimmy SJ. Ren and Li Xu, "On Vectorization of Deep Convolutional Neural Networks for Vision Tasks",
% The 29th AAAI Conference on Artificial Intelligence (AAAI-15). Austin, Texas, USA, January 25-30, 2015
% email: jimmy.sj.ren@gmail.com
%--------------------------------------------------------------------------------------------------------
addpath applications/
addpath utils/
addpath cuda/
addpath mem/
addpath layers/
addpath layers_adapters/
addpath pipeline/
addpath data/MNIST/

clearvars -global config;
clearvars -global mem;

global config;
mnist_configure();
init(0);

%verify_mode = 'gradient';
verify_mode = 'speed';

if strcmp(verify_mode, 'gradient')
    % for gradient checking
    image = config.NEW_MEM(rand(32, 32, 1, config.batch_size)) / (32*32);
    label = config.NEW_MEM((ones(10, config.batch_size) * 2 - 1)) ./ 10;
    op_train_pipe(image, label);
    computeNumericalGradient(image, label, 1);
elseif strcmp(verify_mode, 'speed')
    % for speed testing
    I = gpuArray(single(rand(32,32,1,config.batch_size)));
    y = gpuArray(single(rand(10, config.batch_size)));
    loop = 100;
    tic
    for m = 1:loop
        op_train_pipe(I, y);
    end
    elapse = toc/loop/config.batch_size;
    num = 1 / elapse;
    fprintf('Process %f samples per second.\n', num);
end

author: jimmy-ren


text:
function gptrain_options = get_GPTrainOptions(Ns_gp,hypstruct,optimState,stats,options)
%GETGPTRAINOPTIONS Get options for training GP hyperparameters.

iter = optimState.iter;
if iter > 1; rindex = stats.rindex(iter-1); else; rindex = Inf; end

gptrain_options.OutwarpFun = optimState.gpOutwarpfun;
gptrain_options.Thin = options.GPSampleThin;    % MCMC thinning
gptrain_options.InitMethod = options.GPTrainInitMethod;
gptrain_options.TolOpt = options.GPTolOpt;
gptrain_options.TolOptMCMC = options.GPTolOptMCMC;
gptrain_options.Widths = [];

% Get hyperparameter posterior covariance from previous iters
hypcov = GetHypCov(hypstruct,optimState,stats,options);

% Set up MCMC sampler
switch lower(options.GPHypSampler)
    case {'slicesample'}
        gptrain_options.Sampler = 'slicesample';
        if options.GPSampleWidths > 0 && ~isempty(hypcov)
            widthmult = max(options.GPSampleWidths,rindex);
            hypwidths = sqrt(diag(hypcov)');
            gptrain_options.Widths = max(hypwidths,1e-3)*widthmult;
        else
            gptrain_options.Widths = [];
        end
    case {'npv'}
        gptrain_options.Sampler = 'npv';
    case {'mala'}
        gptrain_options.Sampler = 'mala';
        if ~isempty(hypcov)
            gptrain_options.Widths = sqrt(diag(hypcov)');
        else
            gptrain_options.Widths = [];
        end
        if isfield(optimState,'gpmala_stepsize')
            gptrain_options.Stepsize = optimState.gpmala_stepsize;
        end
    case {'slicelite'}
        gptrain_options.Sampler = 'slicelite';
        if options.GPSampleWidths > 0 && ~isempty(hypcov)
            widthmult = max(options.GPSampleWidths,rindex);
            hypwidths = sqrt(diag(hypcov)');
            gptrain_options.Widths = max(hypwidths,1e-3)*widthmult;
        else
            gptrain_options.Widths = [];
        end
    case {'splitsample'}
        gptrain_options.Sampler = 'splitsample';
        if options.GPSampleWidths > 0 && ~isempty(hypcov)
            widthmult = max(options.GPSampleWidths,rindex);
            hypwidths = sqrt(diag(hypcov)');
            gptrain_options.Widths = max(hypwidths,1e-3)*widthmult;
        else
            gptrain_options.Widths = [];
        end
    case 'covsample'
        if options.GPSampleWidths > 0 && ~isempty(hypcov)
            widthmult = max(options.GPSampleWidths,rindex);
            if all(isfinite(widthmult)) && all(rindex < options.CovSampleThresh)
                nhyp = size(hypcov,1);
                gptrain_options.Widths = (hypcov + 1e-6*eye(nhyp))*widthmult^2;
                gptrain_options.Sampler = 'covsample';
                gptrain_options.Thin = gptrain_options.Thin*ceil(sqrt(nhyp));
            else
                hypwidths = sqrt(diag(hypcov)');
                gptrain_options.Widths = max(hypwidths,1e-3)*widthmult;
                gptrain_options.Sampler = 'slicesample';
            end
        else
            gptrain_options.Widths = [];
            gptrain_options.Sampler = 'slicesample';
        end
    case 'laplace'
        gptrain_options.Widths = [];
        if optimState.Neff < 30
            gptrain_options.Sampler = 'slicesample';
            if options.GPSampleWidths > 0 && ~isempty(hypcov)
                widthmult = max(options.GPSampleWidths,rindex);
                hypwidths = sqrt(diag(hypcov)');
                gptrain_options.Widths = max(hypwidths,1e-3)*widthmult;
            end
        else
            gptrain_options.Sampler = 'laplace';
        end

    otherwise
        error('vbmc:UnknownSampler', ...
            'Unknown MCMC sampler for GP hyperparameters.');
end

% N-dependent initial training points
a = -(options.GPTrainNinit - options.GPTrainNinitFinal);
b = -3*a;
c = 3*a;
d = options.GPTrainNinit;
x = (optimState.Neff - options.FunEvalStart) / (min(options.MaxFunEvals,1e3)-options.FunEvalStart);
f = @(x) a*x.^3 + b*x.^2 + c*x + d;
Ninit = max(round(f(x)),0);

% Set other hyperparameter fitting parameters
if optimState.RecomputeVarPost
    gptrain_options.Burnin = gptrain_options.Thin*Ns_gp;
    gptrain_options.Ninit = Ninit;
    if Ns_gp > 0; gptrain_options.Nopts = 1; else; gptrain_options.Nopts = 2; end
else
    gptrain_options.Burnin = gptrain_options.Thin*3;
    if iter > 1 && stats.rindex(iter-1) < options.GPRetrainThreshold
        gptrain_options.Ninit = 0;
        if strcmpi(options.GPHypSampler,'slicelite')
            gptrain_options.Burnin = max(1,ceil(gptrain_options.Thin*log(stats.rindex(iter-1))/log(options.GPRetrainThreshold)))*Ns_gp;
            gptrain_options.Thin = 1;
        end
        if Ns_gp > 0; gptrain_options.Nopts = 0; else; gptrain_options.Nopts = 1; end
    else
        gptrain_options.Ninit = Ninit;
        if Ns_gp > 0; gptrain_options.Nopts = 1; else; gptrain_options.Nopts = 2; end
    end
end

%gptrain_options.Burnin = 1000;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hypcov = GetHypCov(hypstruct,optimState,stats,options)
%GETHYPCOV Get hyperparameter posterior covariance

if optimState.iter > 1
    if options.WeightedHypCov
        w_list = [];
        hyp_list = [];
        w = 1;
        for i = 1:optimState.iter-1
            if i > 1
                % diff_mult = max(1,log(stats.rindex(optimState.iter-i+1)));
                diff_mult = max(1, ...
                    log(stats.sKL(optimState.iter-i+1) ./ (options.TolsKL*options.FunEvalsPerIter)));
                w = w*(options.HypRunWeight^(options.FunEvalsPerIter*diff_mult));
            end
            if w < options.TolCovWeight; break; end     % Weight is getting too small, break

            hyp = stats.gpHypFull{optimState.iter-i};
            nhyp = size(hyp,2);
            if isempty(hyp_list) || size(hyp_list,2) == size(hyp,1)
                hyp_list = [hyp_list; hyp'];
                w_list = [w_list; w*ones(nhyp,1)/nhyp];
            end
        end

        w_list = w_list / sum(w_list);                  % Normalize weights
        mustar = sum(bsxfun(@times,w_list,hyp_list),1); % Weighted mean

        % Weighted covariance matrix
        nhyp = size(hyp_list,2);
        hypcov = zeros(nhyp,nhyp);
        for j = 1:size(hyp_list,1)
            hypcov = hypcov + ...
                w_list(j)*(hyp_list(j,:)-mustar)'*(hyp_list(j,:)-mustar);
        end
        hypcov = hypcov/(1-sum(w_list.^2));

    else
        hypcov = hypstruct.runcov;
    end
else
    hypcov = [];
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function rez = merge_posthoc2(rez)
%fracse = 0.1;
mu = rez.mu;
fracse = rez.ops.fracse;

ops = rez.ops;
LAM = ops.lam(3) * (20./mu).^2;
Nfilt = rez.ops.Nfilt;

Cmerge = Inf *ones(Nfilt);
tfi = rez.iNeigh;
tf = rez.cProj;
clusterIDs = rez.st3(:,2);
%
nSpikes = size(rez.st3,1);

fmax = zeros(nSpikes,1, 'single');
pairs = {};
for testID = 1:Nfilt
    spikesTest = clusterIDs==testID;
%     tfnew = bsxfun(@plus, tf(spikesTest, :), LAM(tfi(:, testID))'.*mu(tfi(:, testID))');
%     tf(spikesTest, :) = bsxfun(@rdivide, tfnew, sqrt(1+LAM(tfi(:, testID)))');

    pp = tfi(:, testID);
    pp(pp==testID) = [];
    pairs{testID} = pp;
    [~, isame] = min( abs(tfi(:, testID)-testID));
    fmax(spikesTest, 1) = tf(spikesTest, isame);
end

%
inewclust = 0;
clear iMegaC
picked = zeros(Nfilt, 1);
% tic
while 1
    [maxseed, iseed] = max(rez.nbins(1:Nfilt) .* (1-picked), [], 1);
%     [maxseed, iseed] = max(mu(1:Nfilt) .* (1-picked), [], 1);
    if maxseed<500
        break;
    end
    picked(iseed) = 1;
    % 21, 69,
    %
%     iseed = 410;

    run_list = [iseed];
    pair_list = pairs{iseed};
    strun = find(clusterIDs==iseed);

    while ~isempty(pair_list)
        %
%         picked_pairs = rez.nbins(pair_list);

        [mmax, ipair] = max(rez.nbins(pair_list));

        if mmax<100
            break;
        end

        ipair = pair_list(ipair);

        %
        imm = ismember(tfi(:, ipair), run_list);
        if sum(imm)
            %
            new_spikes = find(clusterIDs==ipair);
            f1new = max(tf(new_spikes, imm), [], 2);

            f2new = fmax(new_spikes);

            f1old = fmax(strun);
            f2old = NaN * ones(numel(f1old), 1, 'single');
            i0 = 0;
            for j = 1:length(run_list)
                ifeat = find(tfi(:, run_list(j))==ipair);
                if ~isempty(ifeat)
                    f2old(i0 + (1:rez.nbins(run_list(j))),1) = ...
                        tf(clusterIDs==run_list(j), ifeat);
                    i0 = i0 + rez.nbins(run_list(j));
                end
            end

            f1old(isnan(f2old))=[];
            f2old(isnan(f2old))=[];
            mo = merging_score(f1old - f2old, f1new-f2new, ops.fracse);

            if mo<3
                strun = cat(1, strun, new_spikes);
                run_list(end+1) = ipair;
                picked(ipair)   = 1;
                if mmax>300
                    pair_list = unique(cat(1, pair_list, pairs{ipair}));
                    pair_list(ismember(pair_list, run_list)) = [];
                end
            end
        end
        pair_list(pair_list==ipair) = [];
    end

    inewclust = inewclust + 1;

    iMegaC{inewclust} = run_list;
%     [sum(picked) run_list]
end

% toc
%

iMega = zeros(Nfilt, 1);
for i = 1:length(iMegaC)
   iMega(iMegaC{i}) = iMegaC{i}(1);
end
rez.iMega = iMega;
rez.iMegaC = iMegaC;

rez.st3(:,5) = iMega(rez.st3(:,2));

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
function res = le(a,b)
%LE           Implements  a <= b  for slopes, compares only a.r and b.r
%

% written  12/06/98     S.M. Rump
% modified 09/28/01     S.M. Rump  matrices and multi-dimensional arrays
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if ~isa(a,'slope')
    asize = size(a);
    bsize = b.size;
    res = ( a(:)<=b.r(:,1) );
  elseif ~isa(b,'slope')
    asize = a.size;
    bsize = size(b);
    res = ( a.r(:,1)<=b(:) );
  else
    asize = a.size;
    bsize = b.size;
    res = ( a.r(:,1)<=b.r(:,1) );
  end

  if ~isequal(asize,bsize) & ( prod(asize)~=1 ) & ( prod(bsize)~=1 )
    error('incompatible size for slope <= ')
  end

  if prod(asize)==1
    res = reshape(res,bsize);
  else
    res = reshape(res,asize);
  end

repo: OpenMAS


text:
function fx = p01_fun ( option, nvar, x )

%*****************************************************************************80
%
%% P01_FUN evaluates the function for problem 1.
%
%  Title:
%
%    The Freudenstein-Roth function
%
%  Description:
%
%    One way to use a continuation code as a nonlinear root finder
%    is to start with a set of nonlinear equations G(X), and an
%    approximate root A, and create a "homotopy" function F(X,Y)
%    with the properties that F(A,0.0) = 0 and F(X,1.0) = G(X).
%    Thus, the homotopy function F has a known exact solution
%    from which we can start with no difficulty.  If the continuation
%    code can take us from Y = 0 to Y = 1, then we have found
%    an X so that F(X,1.0) = 0, so we have found a solution to G(X)=0.
%
%    The Freudenstein-Roth function F(X) is derived in this way
%    from a homotopy of G(X):
%
%      F ( X(1), X(2), X(3) ) =
%        G ( X(1), X(2) ) - ( 1 - X(3) ) * G ( Y1, Y2 )
%
%    where Y1 and Y2 are some fixed values, and
%
%      G(1) = X(1) - X(2)*X(2)*X(2) + 5*X(2)*X(2) -  2*X(2) - 13
%      G(2) = X(1) + X(2)*X(2)*X(2) +   X(2)*X(2) - 14*X(2) - 29
%
%  Options 1, 2, 3:
%
%    The starting point is X0 = ( 15, -2, 0 ).
%
%    A great deal of information is available about the homotopy curve
%    generated by this starting point:
%
%    The function F(X) has the form
%
%      F(1) = X(1) - X(2)**3 + 5*X(2)**2 -  2*X(2) - 13 + 34*(X(3)-1)
%      F(2) = X(1) + X(2)**3 +   X(2)**2 - 14*X(2) - 29 + 10*(X(3)-1)
%
%    There is a closed form representation of the curve in terms of the
%    second parameter:
%
%      X(1) = (-11*X(2)**3 + 4*X(2)**2 + 114*X(2) + 214) /  6
%      X(2) = X(2)
%      X(3) = (    X(2)**3 - 2*X(2)**2 -   6*X(2) +   4) / 12
%
%    The first option simply requests the production of solution points
%    along the curve until a point is reached whose third component is
%    exactly 1.
%
%    Options 2 and 3 use the same starting point, and also stop when the
%    third component is 1.  However, these options in addition search
%    for limit points in the first and third components of the solution,
%    respectively.
%
%    The target solution has X(3) = 1, and is ( 5, 4, 1 ).
%
%    Limit points for X1:
%
%      ( 14.28309, -1.741377,  0.2585779 )
%      ( 61.66936,  1.983801, -0.6638797 )
%
%    Limit points for X3:
%
%     (20.48586, -0.8968053, 0.5875873)
%     (61.02031,  2.230139, -0.6863528)
%
%    The curve has several dramatic bends.
%
%
%  Options 4, 5, and 6:
%
%    The starting point is (4, 3, 0).
%
%    The function F(X) has the form
%
%      F(1) = X(1) - X(2)**3 + 5*X(2)**2 -  2*X(2) - 13 +  3*(X(3)-1)
%      F(2) = X(1) + X(2)**3 +   X(2)**2 - 14*X(2) - 29 - 31*(X(3)-1)
%
%    There is a closed form representation of the curve in terms of the
%    second parameter:
%
%      X(1) = (14*X(2)**3 -79*X(2)**2 +52*X(2) + 245) / 17
%      X(2) = X(2)
%      X(3) = (   X(2)**3 - 2*X(2)**2 - 6*X(2) +   9) / 17
%
%    The correct value of the solution at X(3)=1 is:
%
%      (5, 4, 1)
%
%    In option 5, limit points in the first component are sought,
%    and in option 6, limit points in the third component are
%    sought.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ferdinand Freudenstein, Bernhard Roth,
%    Numerical Solutions of Nonlinear Equations,
%    Journal of the Association for Computing Machinery,
%    Volume 10, 1963, Pages 550-556.
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the function.
%
%    Output, real FX(NVAR-1), the value of the function at X.
%

%  Get the starting point, Y.
%
  y = p01_start ( option, nvar );
%
%  G is the function value at the starting point,
%  F the function value at the current point.
%
  gy = p01_gx ( y );

  gx = p01_gx ( x );
%
%  The parameter X3 generates the homotopy curve.
%
  fx(1:nvar-1) = gx(1:nvar-1) + ( x(3) - 1.0 ) * gy(1:nvar-1);

  return
end


text:
%{
This code demonstrates how to warp a depth map (togeter with a label map) from a camera to another view.
It create a mesh using the depth map from a camera and render the mesh in OpenGL.
Therefore, it will not create artifacts with a lot of gaps in a warpped result.

compile
   mex WarpMesh.cpp -lGLU -lOSMesa
   you need to install osmesa to compile it
   in Mac OS X 10.9.2, you just need to install X11 on mac (xquartz) comes with mesa
   mex WarpMesh.cpp -lGLU -lOSMesa -I/opt/X11/include/ -L/opt/X11/lib/

Email Jianxiong if you have questions. Cite the following paper if you use this code:

J. Xiao, A. Owens and A. Torralba
SUN3D: A Database of Big Spaces Reconstructed using SfM and Object Labels
Proceedings of 14th IEEE International Conference on Computer Vision (ICCV2013)

%}

close all;
clear

load demo.mat

[label,depth] = WarpMeshMatlab(XYZcamera,labelNow,K, Rt);

% warp color
[imageWarp,depth] = WarpMeshMatlabColor(XYZcamera,double(image),K, Rt);

figure
subplot(2,3,1);
imagesc(XYZcamera(:,:,3));
axis equal; axis tight; axis off; title('before warping');

subplot(2,3,2);
imagesc(labelNow);
axis equal; axis tight; axis off; title('before warping');

subplot(2,3,4);
imagesc(depth);
axis equal; axis tight; axis off; title('after warping');

subplot(2,3,5);
imagesc(label);
axis equal; axis tight; axis off; title('after warping');

subplot(2,3,3);
imagesc(image);
axis equal; axis tight; axis off; title('before warping');

subplot(2,3,6);
imagesc(imageWarp/255);
axis equal; axis tight; axis off; title('after warping');

author: jianxiongxiao

repo: ProfXkit


text:
function demnlab(action);
%DEMNLAB A front-end Graphical User Interface to the demos
%
%	Description
%	This function will start a user interface allowing the user to select
%	different demonstration functions to view. The demos are divided into
%	4 groups, with the demo being executed by selecting the desired
%	option from a pop-up menu.
%
%	See also
%

%	Copyright (c) Ian T Nabney (1996-2001)

% If run without parameters, initialise gui.
if nargin<1,
	action='initialise';
end;

if strcmp(action,'initialise'),

  % Create figure
  fig = figure( ...
    'Name', 'Netlab Demos', ...
    'NumberTitle', 'off', ...
    'Color', [0.7529 0.7529 0.7529], ...
    'Visible', 'on');

  % Create GROUPS
  % Bottom of demo buttons
  group1_bot = 0.20;
  group1_top = 0.75;
  uicontrol(fig, ...
    'Style', 'frame', ...
    'Units', 'normalized', ...
    'Position', [0.03 group1_bot 0.94 group1_top - group1_bot], ...
    'BackgroundColor', [0.5 0.5 0.5]);

  % Bottom of help and close buttons
  group2_bot = 0.04;
  uicontrol(fig, ...
    'Style', 'frame', ...
    'Units', 'normalized', ...
    'Position', [0.03 group2_bot 0.94 0.12], ...
    'BackgroundColor', [0.5 0.5 0.5]);

  % Draw title
  hLogoAxis = axes( ...
    'Units', 'normalized', ...
    'Position', [0.05 0.82 0.90 0.14], ...
    'Box', 'off', ...
    'XColor', [0 0 0], ...
    'YColor', [0 0 0], ...
    'Visible', 'on');

  load netlogo;			% load image and colour map
  colormap(netcmap(1:3,:));	% change colour map: don't need many entries
  image(nlogo);			% draw logo
  axis('image');		% ensures pixels on axis are square
  axis off;			% turn axes off

  % Create static text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.5 0.5 0.5], ...
    'Position', [0.05 group1_top-0.1 0.90 0.08], ...
    'String', 'Select demo to run:');

  % First row text offset
  tRow1Offset = 0.14;
  % Offset between text and button
  TBoffset = 0.07;
  % First row button offset
  bRow1Offset = tRow1Offset+TBoffset;
  % ONE text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.5 0.5 0.5], ...
    'Position', [0.08 group1_top-tRow1Offset 0.36 0.05], ...
    'String', 'Regression');

  popup1str(1) = {'Select Option'};
  popup1str(2) = {'Multi-Layer Perceptron'};
  popup1str(3) = {'Radial Basis Function'};
  popup1str(4) = {'Mixture Density Network'};
  % ONE popup
  hPop1 = uicontrol(fig, ...
    'Style','popup', ...
    'Units','normalized', ...
    'String', popup1str, ...
    'Position', [0.08 group1_top-bRow1Offset 0.36 0.08], ...
    'BackgroundColor', [0.7 0.7 0.7], ...
    'Callback', 'demnlab popup1');

  % TWO text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.5 0.5 0.5], ...
    'Position', [0.56 group1_top-tRow1Offset 0.36 0.05], ...
    'String', 'Classification');

  popup2str(1) = popup1str(1);
  popup2str(2) = {'Generalised Linear Model (2 class)'};
  popup2str(3) = {'Generalised Linear Model (3 class)'};
  popup2str(4) = {'Multi-Layer Perceptron'};
  popup2str(5) = {'K nearest neighbour'};
  % TWO popup
  hPop2 = uicontrol(fig, ...
    'Style','popup', ...
    'Units','normalized', ...
    'String', popup2str, ...
    'Position', [0.56 group1_top-bRow1Offset 0.36 0.08], ...
    'BackgroundColor', [0.7 0.7 0.7], ...
    'Callback', 'demnlab popup2');

  tRow2Offset = 0.30;
  bRow2Offset = tRow2Offset+TBoffset;
  % THREE text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.5 0.5 0.5], ...
    'Position', [0.08 group1_top - tRow2Offset 0.36 0.05], ...
    'String', 'Density Modelling and Clustering');

  popup3str(1) = popup1str(1);
  popup3str(2) = {'Gaussian Mixture (EM training)'};
  popup3str(3) = {'Gaussian Mixture (spherical)'};
  popup3str(4) = {'Gaussian Mixture (diagonal)'};
  popup3str(5) = {'Gaussian Mixture (full)'};
  popup3str(6) = {'Neuroscale'};
  popup3str(7) = {'GTM (EM training)'};
  popup3str(8) = {'GTM (visualisation)'};
  popup3str(9) = {'K-means clustering'};
  popup3str(10) = {'Self-Organising Map'};
  % TWO popup
  % THREE popup
  hPop3 = uicontrol(fig, ...
    'Style','popup', ...
    'Units','normalized', ...
    'String', popup3str, ...
    'Position', [0.08 group1_top - bRow2Offset 0.36 0.08], ...
    'BackgroundColor', [0.7 0.7 0.7], ...
    'Callback', 'demnlab popup3');

  % FOUR text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.5 0.5 0.5], ...
    'Position', [0.56 group1_top - tRow2Offset 0.36 0.05], ...
    'String', 'Bayesian Methods');

  popup4str(1) = popup1str(1);
  popup4str(2) = {'Sampling the MLP Prior'};
  popup4str(3) = {'Evidence Approximation for MLP'};
  popup4str(4) = {'Evidence Approximation for RBF'};
  popup4str(5) = {'Evidence Approximation in Classification'};
  popup4str(6) = {'ARD for MLP'};
  popup4str(7) = {'Sampling the GP Prior'};
  popup4str(8) = {'GPs for Regression'};
  popup4str(9) = {'ARD for GP'};
  % FOUR popup
  hPop4 = uicontrol(fig, ...
    'Style','popup', ...
    'Units','normalized', ...
    'String', popup4str, ...
    'Position', [0.56 group1_top - bRow2Offset 0.36 0.08], ...
    'BackgroundColor', [0.7 0.7 0.7], ...
    'Callback', 'demnlab popup4');

  tRow3Offset = 0.45;
  bRow3Offset = tRow3Offset+TBoffset;
  % FIVE text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.5 0.5 0.5], ...
    'Position', [0.08 group1_top - tRow3Offset 0.36 0.05], ...
    'String', 'Optimisation and Visualisation');

  popup5str(1) = popup1str(1);
  popup5str(2) = {'Algorithm Comparison'};
  popup5str(3) = {'On-line Gradient Descent'};
  popup5str(4) = {'Hinton Diagrams'};
  % FIVE popup
  hPop5 = uicontrol(fig, ...
    'Style','popup', ...
    'Units','normalized', ...
    'String',popup5str, ...
    'Position', [0.08 group1_top - bRow3Offset 0.36 0.08], ...
    'BackgroundColor', [0.7 0.7 0.7], ...
    'Callback', 'demnlab popup5');

  % SIX text
  uicontrol(fig, ...
    'Style', 'text', ...
    'Units', 'normalized', ...
    'BackgroundColor', [0.5 0.5 0.5], ...
    'Position', [0.56 group1_top - tRow3Offset 0.36 0.05], ...
    'String', 'Sampling');

  popup6str(1) = popup1str(1);
  popup6str(2) = {'Sampling a Gaussian'};
  popup6str(3) = {'MCMC sampling (Metropolis)'};
  popup6str(4) = {'Hybrid MC (Gaussian mixture)'};
  popup6str(5) = {'Hybrid MC for MLP I'};
  popup6str(6) = {'Hybrid MC for MLP II'};
  % SIX popup
  hPop6 = uicontrol(fig, ...
    'Style','popup', ...
    'Units','normalized', ...
    'String', popup6str, ...
    'Position', [0.56 group1_top - bRow3Offset 0.36 0.08], ...
    'BackgroundColor', [0.7 0.7 0.7], ...
    'Callback', 'demnlab popup6');

  % Create HELP button
  uicontrol(fig, ...
    'Units', 'normalized', ...
    'Position' , [0.05 group2_bot+0.02 0.40 0.08], ...
    'String', 'Help', ...
    'Callback', 'demnlab help');

  % Create CLOSE button
  uicontrol(fig, ...
    'Units', 'normalized', ...
    'Position' , [0.55 group2_bot+0.02 0.40 0.08], ...
    'String', 'Close', ...
    'Callback', 'close(gcf)');

  hndlList=[fig hPop1 hPop2 hPop3 hPop4 hPop5 hPop6];
  set(fig, 'UserData', hndlList);
  set(fig, 'HandleVisibility', 'callback');

elseif strcmp(action, 'popup1'),

  hndlList=get(gcf,'UserData');
  hPop = hndlList(2);

  selected = get(hPop, 'Val');
  set(hPop, 'Val', [1]);

  switch selected
    case 2
      demmlp1;
    case 3
      demrbf1;
    case 4
      demmdn1;
  end;

elseif strcmp(action,'popup2'),

  hndlList=get(gcf,'UserData');
  hPop = hndlList(3);

  selected = get(hPop, 'Val');
  set(hPop, 'Val', [1]);

  switch selected
    case 2
      demglm1;
    case 3
      demglm2;
    case 4
      demmlp2;
    case 5
      demknn1;
  end

elseif strcmp(action,'popup3'),

  hndlList=get(gcf,'UserData');
  hPop = hndlList(4);

  selected = get(hPop, 'Val');
  set(hPop, 'Val', [1]);

  switch selected
    case 2
      demgmm1;
    case 3
      demgmm2;
    case 4
      demgmm3;
    case 5
      demgmm4;
    case 6
      demns1;
    case 7
      demgtm1;
    case 8
      demgtm2;
    case 9
      demkmn1;
    case 10
      demsom1;
  end

elseif strcmp(action,'popup4'),

  hndlList=get(gcf,'UserData');
  hPop = hndlList(5);

  selected = get(hPop, 'Val');
  set(hPop, 'Val', [1]);

  switch selected
    case 2
      demprior;
    case 3
      demev1;
  case 4
      demev3;
  case 5
      demev2;
    case 6
      demard;
    case 7
      demprgp;
    case 8
      demgp;
    case 9
      demgpard;
  end

elseif strcmp(action,'popup5'),

  hndlList=get(gcf,'UserData');
  hPop = hndlList(6);

  selected = get(hPop, 'Val');
  set(hPop, 'Val', [1]);

  switch selected
    case 2
      demopt1;
    case 3
      demolgd1;
    case 4
      demhint;
  end

elseif strcmp(action,'popup6'),

  hndlList=get(gcf,'UserData');
  hPop = hndlList(7);

  selected = get(hPop, 'Val');
  set(hPop, 'Val', [1]);

  switch selected
    case 2
      demgauss;
    case 3
      demmet1;
    case 4
      demhmc1;
    case 5
      demhmc2;
    case 6
      demhmc3;
  end

elseif strcmp(action, 'help'),

  helpStr = {'To run a demo, press the appropriate button.'; ...
	'Instructions and information will appear in the Matlab';...
	'command window.'};

  hHelpDlg = helpdlg(helpStr, 'Netlab Demo Help');

end;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% corrcoef_cell() - compute pairwise correlations using arrays and
%                   cell array inputs.
%
% Usage:
%    >> c = corrcoef_cell( data );
%    >> c = corrcoef_cell( data );
%
% Inputs:
%   data       - [cell array] data consisting of PAIRED arrays to be compared.
%                The last dimension of embeded data arrays is used to compute
%                correlation (see examples).
% Outputs:
%   c   - Correlation values. Same size as data without the last dimension.
%
% Note: the main advantage over the corrcoef Matlab function is the
%       capacity to compute millions of pairwise correlations per second.
%
% Example:
%   a = { rand(1,10) rand(1,10) };
%   c1 = corrcoef_cell(a);
%   c2 = corrcoef(a{1}, a{2});
%   % in this case, c1 is equal to c2(2)
%
%   a = { rand(200,300,100) rand(200,300,100) };
%   c = corrcoef_cell(a);
%   % the call above would require 200 x 300 calls to the corrcoef function
%   % and be about 1000 times slower
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2010

% Copyright (C) Arnaud Delorme
%
% 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 c = corrcoef_cell(a,b);

if nargin < 1
    help corrcoef_cell;
    return;
end;

if nargin < 2
    b = a{2};
    a = a{1};
end;

nd = myndims(a);
if nd == 1
    aa = a-mean(a);
    bb = b-mean(b);
    cv  = aa'*bb/(10-1);
    cva = aa'*aa/(10-1);
    cvb = bb'*bb/(10-1);

    c = cv/sqrt(cva*cvb);
elseif nd == 2
    aa = a-repmat(mean(a,2),[1 size(a,2)]);
    bb = b-repmat(mean(b,2),[1 size(a,2)]);
    cv  = sum(aa.*bb,2);
    cva = sum(aa.*aa,2);
    cvb = sum(bb.*bb,2);

    c = cv./sqrt(cva.*cvb);
elseif nd == 3
    aa = a-repmat(mean(a,3),[1 1 size(a,3)]);
    bb = b-repmat(mean(b,3),[1 1 size(a,3)]);
    cv  = sum(aa.*bb,3);
    cva = sum(aa.*aa,3);
    cvb = sum(bb.*bb,3);

    c = cv./sqrt(cva.*cvb);
elseif nd == 4
    aa = a-repmat(mean(a,4),[1 1 1 size(a,4)]);
    bb = b-repmat(mean(b,4),[1 1 1 size(a,4)]);
    cv  = sum(aa.*bb,4);
    cva = sum(aa.*aa,4);
    cvb = sum(bb.*bb,4);

    c = cv./sqrt(cva.*cvb);
end;

function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end;
    end;


text:
%
% [K, dK] = gpK_decper(D, p1, p2, p3, p4)
%
% Decaying periodic covariance function for GP.
%
% p1 is the scale
% p2 is the decaying length
% p3 is the frequency of the periodicity
% p4 is the smoothness of the periodic signal
%
function [K, dK] = gpK_decper(D, p1, p2, p3, p4)
K = p1^2*exp(-0.5*(D.^2)/(p2^2)-2*sin(pi*D*p3).^2/(p4^2));

if nargout >= 2
  dK = zeros([size(D),4]);
  dK(:,:,1) = K .* 2 / p1;
  dK(:,:,2) = K .* (-0.5*D.^2) .* (-2*p2^(-3));
  dK(:,:,3) = K .* (-2*sin(pi*D*p3)*2) .* cos(pi*D*p3) .* (pi*D);
  dK(:,:,4) = K .* (-2*sin(pi*D*p3).^2) * (-2)*p4^(-3);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [ ret_box8d ] = getRotateBox8D( box8d, orientation, center_x, center_y)
    boxes = box8d;
    rotate_mat = [cosd(orientation), -sind(orientation); sind(orientation), cosd(orientation)];

    %% move to center
    boxes(:,1) = boxes(:,1) - center_y;
    boxes(:,2) = boxes(:,2) - center_x;
    boxes = boxes(:,[2 1])';
    ret_box8d = rotate_mat * boxes;
    ret_box8d = ret_box8d([2,1],:)';

    %% resotre position
    ret_box8d(:,1) = ret_box8d(:,1) + center_y;
    ret_box8d(:,2) = ret_box8d(:,2) + center_x;

    ret_box8d = round(ret_box8d);
end

author: stupidZZ


text:
%% Simulate some RF pulses and analyze their properties
%
% This example demonstrates the use of mr.simRf() function.

%% Create a system object
% we need it here to create fat-sat pulses
sys = mr.opts('B0', 2.89);
%seq=mr.Sequence(sys);

%% 30 degree slice selective SINC pulse
rf30_sinc = mr.makeSincPulse(pi/6,'system',sys,'Duration',3e-3,'use','excitation',...
    'PhaseOffset',pi/2,'apodization',0.3,'timeBwProduct',4);

[bw,f0,M_xy_sta,F1]=mr.calcRfBandwidth(rf30_sinc);
[M_z,M_xy,F2]=mr.simRf(rf30_sinc);

figure; plot(F1,abs(M_xy_sta),F2,abs(M_xy),F2,M_z);
axis([f0-2*bw, f0+2*bw, -0.1, 1.2]);
legend({'M_x_ySTA','M_x_ySIM','M_zSIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('STA vs. simulation, flip angle 30°');

figure; plot(F2,real(M_xy),F2,imag(M_xy));
axis([f0-2*bw, f0+2*bw, -1.2, 1.2]);
legend({'M_xSIM','M_ySIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('Real and imag. parts of transverse magnetisation, 30° flip');

%% 90 degree slice selective SINC pulse
rf90_sinc = mr.makeSincPulse(pi/2,'system',sys,'Duration',3e-3,'use','excitation',...
    'PhaseOffset',pi/2,'apodization',0.6,'timeBwProduct',8);

[bw,f0,M_xy_sta,F1]=mr.calcRfBandwidth(rf90_sinc);
[M_z,M_xy,F2]=mr.simRf(rf90_sinc);

%%
figure; plot(F1,abs(M_xy_sta),F2,abs(M_xy),F2,M_z);
axis([f0-2*bw, f0+2*bw, -0.1, 1.2]);
legend({'M_x_ySTA','M_x_ySIM','M_zSIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('STA vs. simulation, flip angle 90°');

figure; plot(F2,atan2(abs(M_xy),M_z)/pi*180);
axis([f0-2*bw, f0+2*bw, -5, 100]);
xlabel('frequency offset / Hz');
ylabel('flip ange [°]');
legend({'SINC'});
grid on;
title('Achieved flip angle for the nominal 90° flip');

figure; plot(F2,real(M_xy),F2,imag(M_xy));
axis([f0-2*bw, f0+2*bw, -1.2, 1.2]);
legend({'M_xSIM','M_ySIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('Real and imag. parts of transverse magnetisation, 90° flip');

%%
figure; plot(F2,angle(M_xy));
axis([f0-2*bw, f0+2*bw, -3.2, 3.2]);
xlabel('frequency offset / Hz');
ylabel('phase');
title('Phase of transverse magnetisation, 90° flip');

M_xy_masked=M_xy;
M_xy_masked(F2<f0-bw/2)=0;
M_xy_masked(F2>f0+bw/2)=0;

i_phase_slope=sum(M_xy_masked(2:end).*conj(M_xy_masked(1:end-1)));
i_phase_slope=1i*angle(i_phase_slope)/(F2(2)-F2(1));
i_phase_offset=sum(M_xy_masked.*exp(-i_phase_slope*F2));
i_phase_offset=i_phase_offset/abs(i_phase_offset);

hold on; plot(F2,angle(exp(i_phase_slope*F2)*i_phase_offset),'--');

xline(f0,'-');
xline(f0-bw/2,'--');
xline(f0+bw/2,'--');

legend({'SINC-phase','linear fit'});
fprintf('SINC90 rf center error: %g (%g %%)\n', abs(i_phase_slope)/2*pi/rf90_sinc.shape_dur/10, 100/0.5*abs(i_phase_slope)/2*pi/rf90_sinc.shape_dur/10); % no idea where this 10 comes from

%% 90 degree slice selective SLR pulse
rf_90slr= mr.makeSLRpulse(pi/2,'duration',3e-3,'timeBwProduct',4,'PhaseOffset',pi/2,'use','excitation',...
    'passbandRipple',1,'stopbandRipple',1e-2,'filterType','ms','system',sys);

[bw,f0,M_xy_sta,F1]=mr.calcRfBandwidth(rf_90slr);
[M_z,M_xy,F2]=mr.simRf(rf_90slr);

figure; plot(F1,abs(M_xy_sta),F2,abs(M_xy),F2,M_z);
axis([f0-2*bw, f0+2*bw, -0.1, 1.2]);
legend({'M_x_ySTA','M_x_ySIM','M_zSIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('SLR: STA vs. simulation, flip angle 90°');

figure; plot(F2,real(M_xy),F2,imag(M_xy));
axis([f0-2*bw, f0+2*bw, -1.2, 1.2]);
legend({'M_xSIM','M_ySIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('SLR: Real and imag. parts of transverse magnetisation, 90° flip');

figure; plot(F2,angle(M_xy));
axis([f0-2*bw, f0+2*bw, -3.2, 3.2]);
xlabel('frequency offset / Hz');
ylabel('phase');
legend({'SLR'});
title('SLR: Phase of transverse magnetisation, 90° flip');

figure; plot(F2,atan2(abs(M_xy),M_z)/pi*180);
axis([f0-2*bw, f0+2*bw, -5, 100]);
xlabel('frequency offset / Hz');
ylabel('flip ange [°]');
legend({'SLR'});
grid on;
title('SLR: Achieved flip angle for the nominal 90° flip');

%% 60 degree slice selective SINC pulse
rf60_block = mr.makeBlockPulse(pi/3,'system',sys,'Duration',0.5e-3,'use','excitation','PhaseOffset',pi/2);

[bw,f0,M_xy_sta,F1]=mr.calcRfBandwidth(rf60_block);
[M_z,M_xy,F2]=mr.simRf(rf60_block);

figure; plot(F1,abs(M_xy_sta),F2,abs(M_xy),F2,M_z);
axis([f0-2*bw, f0+2*bw, -0.1, 1.2]);
legend({'M_x_ySTA','M_x_ySIM','M_zSIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('STA vs. simulation, 60° hard pulse');

figure; plot(F2,real(M_xy),F2,imag(M_xy));
axis([f0-2*bw, f0+2*bw, -1.2, 1.2]);
legend({'M_xSIM','M_ySIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('Real and imag. parts of transverse magnetisation, 60° hard pulse');

%% fat-sat pulse
sat_ppm=-3.45;
sat_freq=sat_ppm*1e-6*sys.B0*sys.gamma;
rf_fs = mr.makeGaussPulse(110*pi/180,'system',sys,'Duration',8e-3,...
    'bandwidth',abs(sat_freq),'freqOffset',sat_freq,'use','saturation');
rf_fs.phaseOffset=-2*pi*rf_fs.freqOffset*mr.calcRfCenter(rf_fs); % compensate for the frequency-offset induced phase

[M_z,M_xy,F2]=mr.simRf(rf_fs);

figure; plot(F2,real(M_xy),F2,imag(M_xy),F2,M_z);
axis([sat_freq-900, sat_freq+900, -1.2, 1.2]);
legend({'M_x','M_y','M_z'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('Simulation, Gaussian fat-sat pulse');

%% 180 degree slice selective SINC pulse
rf180_sinc = mr.makeSincPulse(pi,'system',sys,'Duration',4e-3,'use','refocusing',...
    'apodization',0.3,'timeBwProduct',6);

[bw,f0,M_xy_sta,F1]=mr.calcRfBandwidth(rf180_sinc);
[M_z,M_xy,F2,ref_eff]=mr.simRf(rf180_sinc);

% figure; plot(F1,abs(M_xy_sta),F2,abs(M_xy),F2,M_z);
% axis([f0-2*bw, f0+2*bw, -1.2, 1.2]);
% legend({'M_x_ySTA','M_x_ySIM','M_zSIM'});
% xlabel('frequency offset / Hz');
% ylabel('magnetisation');
% title('STA vs. simulation, flip angle 180°');
%
% figure; plot(F2,real(M_xy),F2,imag(M_xy));
% axis([f0-2*bw, f0+2*bw, -1.2, 1.2]);
% legend({'M_xSIM','M_ySIM'});
% xlabel('frequency offset / Hz');
% ylabel('magnetisation');
% title('Real and imag. parts of transverse magnetisation, 180° flip');
%
% figure; plot(F2,angle(M_xy));
% axis([f0-2*bw, f0+2*bw, -3.2, 3.2]);
% xlabel('frequency offset / Hz');
% ylabel('phase');
% title('Phase of transverse magnetisation, 180° flip');

figure; plot(F2,atan2(abs(M_xy),M_z)/pi*180);
axis([f0-2*bw, f0+2*bw, -5, 190]);
xlabel('frequency offset / Hz');
ylabel('flip ange [°]');
legend({'SINC'});
grid on;
title('Achieved flip angle for the nominal 180° flip');

figure; plot(F2,abs(ref_eff));
axis([f0-2*bw, f0+2*bw, -0.1, 1.1]);
xlabel('frequency offset / Hz');
ylabel('efficiency');
legend({'SINC'});
title('refocusing efficiency');

figure; plot(F2,angle(ref_eff));
axis([f0-2*bw, f0+2*bw, -3.2, 3.2]);
xlabel('frequency offset / Hz');
ylabel('efficiency');
legend({'SINC'});
title('refocusing efficiency phase (~2x RF phase)');

%% 180 degree slice selective SLR pulse
rf180_slr= mr.makeSLRpulse(pi,'duration',4e-3,'timeBwProduct',6,'use','refocusing','filterType','ms','system',sys);

[bw,f0,M_xy_sta,F1]=mr.calcRfBandwidth(rf180_slr);
[M_z,M_xy,F2_slr,ref_eff_slr]=mr.simRf(rf180_slr);

figure; plot(F2_slr,atan2(abs(M_xy),M_z)/pi*180);
axis([f0-2*bw, f0+2*bw, -5, 190]);
xlabel('frequency offset / Hz');
ylabel('flip ange [°]');
legend({'SLR'});
grid on;
title('Achieved flip angle for the nominal 180° flip');

figure; plot(F2_slr,abs(ref_eff_slr));
axis([f0-2*bw, f0+2*bw, -0.1, 1.1]);
xlabel('frequency offset / Hz');
ylabel('efficiency');
legend({'SLR'});
title('refocusing efficiency');

figure; plot(F2_slr,angle(ref_eff_slr));
axis([f0-2*bw, f0+2*bw, -3.2, 3.2]);
xlabel('frequency offset / Hz');
ylabel('efficiency');
title('SLR: refocusing efficiency phase (~2x RF phase)');

figure; plot(F2,abs(ref_eff),F2_slr,abs(ref_eff_slr));
axis([f0-2*bw, f0+2*bw, -0.1, 1.1]);
xlabel('frequency offset / Hz');
ylabel('efficiency');
legend({'SINC','SLR'});
title('refocusing efficiency: SINC vs SLR');

%% adiabatic pulse
rf180_ad = mr.makeAdiabaticPulse('wurst','duration',4.6e-3,'bandwidth',6000,'n_fac',20,'use','inversion','system',sys);

[bw,f0,M_xy_sta,F1]=mr.calcRfBandwidth(rf180_ad);
[M_z,M_xy,F2,ref_eff,ref_mx,ref_my]=mr.simRf(rf180_ad,-0.5);

figure; plot(F1,abs(M_xy_sta),F2,abs(M_xy),F2,M_z);
axis([f0-2*bw, f0+2*bw, -1.2, 1.2]);
legend({'M_x_ySTA','M_x_ySIM','M_zSIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('STA vs. simulation, adiabatic pulse');

figure; plot(F2,real(M_xy),F2,imag(M_xy));
axis([f0-2*bw, f0+2*bw, -1.2, 1.2]);
legend({'M_xSIM','M_ySIM'});
xlabel('frequency offset / Hz');
ylabel('magnetisation');
title('Real and imag. parts of transverse magnetisation, adiabatic pulse');

figure; plot(F2,angle(M_xy));
axis([f0-2*bw, f0+2*bw, -3.2, 3.2]);
xlabel('frequency offset / Hz');
ylabel('phase');
legend({'WURST'});
title('Phase of transverse magnetisation, adiabatic pulse');

figure; plot(F2,atan2(abs(M_xy),M_z)/pi*180);
axis([f0-2*bw, f0+2*bw, -5, 190]);
xlabel('frequency offset / Hz');
ylabel('flip ange [°]');
legend({'WURST'});
grid on;
title('Achieved flip angle adiabatic pulse');

figure; plot(F2,abs(ref_eff));
axis([f0-2*bw, f0+2*bw, -0.1, 1.1]);
xlabel('frequency offset / Hz');
ylabel('efficiency');
legend({'WURST'});
title('refocusing efficiency');

figure; plot(F2,angle(ref_eff));
axis([f0-2*bw, f0+2*bw, -3.2, 3.2]);
xlabel('frequency offset / Hz');
ylabel('efficiency');
legend({'WURST'});
title('refocusing efficiency phase (~2x RF phase)');

%% spoiling simulation for the same pulse used and refocusing pulse

spoiling_factor=5; % area of the left/righ spoiler; reasonable range 1..10
cl=13; % convolution length to simulate intravoxel dephasing

[M_z,M_xy,F2,ref_eff,mxrf,myrf]=mr.simRf(rf180_ad,spoiling_factor,spoiling_factor);

mxrfc=conv(mxrf,ones(cl,1)/cl,'same');
myrfc=conv(myrf,ones(cl,1)/cl,'same');

figure;plot(F2,abs(abs(mxrfc)+1i*abs(myrfc))/2^0.5,F2,abs(ref_eff),F2,0.5-0.5*M_z,'--');
xlabel('frequency offset / Hz');
ylabel('signal');
legend({'spoiling','ref.eff.','.5-.5*M_z'});
title('signal with spoiling vs refocusing efficiency');

%% investigare RF center shift as a function of clip angle (is it an artifact?)

alphas=[5:5:150];
rfce=[];
for a=alphas
    rfAlpha = mr.makeSincPulse(a/180*pi,'system',sys,'Duration',3e-3,'use','excitation','apodization',0.3,'timeBwProduct',4);
    %rfAlpha = mr.makeGaussPulse(a*pi/180,'system',sys,'Duration',3e-3,'timeBwProduct',8,'use','excitation');

    [bw,f0,M_xy_sta,F1]=mr.calcRfBandwidth(rfAlpha);
    [M_z,M_xy,F2]=mr.simRf(rfAlpha);

    M_xy_masked=M_xy;
    M_xy_masked(F2<f0-bw/2)=0;
    M_xy_masked(F2>f0+bw/2)=0;

    i_phase_slope=sum(M_xy_masked(2:end).*conj(M_xy_masked(1:end-1)));
    i_phase_slope=1i*angle(i_phase_slope)/(F2(2)-F2(1));

    rfce=[rfce,abs(i_phase_slope)/2*pi/rf90_sinc.shape_dur/10]; % no idea where this 10 comes from
end

figure;plot(alphas,rfce/0.5*100); title('excess gradient refocusing needed in %');
xticks([0:30:alphas(end)]);
yticks([0:2:10]);
grid on;
xlabel('flip angle / °');
ylabel('refocusing / %');

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function [b, r, crit] = stepwise(X, y, criterion)

    if nargin<3, criterion = 'BIC'; end

    % qr factorization will speed up stepwise regression significantly
    [Q,R] = qr(X,0); % note that the zero is very important for performance
    invR = pinv(R);

    n = size(y,1);
    LL = nan(size(X,2),1);
    for i = 1:length(LL)
        % get residual for each fit
        b = invR(1:i,1:i) * Q(:,1:i)' * y;
        r = y - X(:,1:i)*b;

        % calculate log-likelihood
        LL(i) = -n/2*log( 2*pi*mean(r.^2) ) - n/2;

    end

    % Calculate information criterion
    crit = infocrit( LL , n , (1:length(LL))' , criterion );

    % optimal model order
    lst=find(~isnan(crit));
    [~, N] = min( crit(lst) );
    N=lst(N);

    % finally, our output
    b = invR(1:N,1:N) * Q(:,1:N)'*y;
    r = y - X(:,1:N)*b;

end

author: BUNPC


text:
function variance = nakagami_variance ( a, b, c )

%*****************************************************************************80
%
%% NAKAGAMI_VARIANCE returns the variance of the Nakagami PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < B
%    0.0 < C
%
%    Output, real VARIANCE, the variance of the PDF.
%
  t1 = gamma ( c + 0.5 );
  t2 = gamma ( c );

  variance = b * b * ( 1.0 - t1 * t1 / ( c * t2 * t2 ) );

  return
end


text:
function [val] = optGeneFitnessTilt(rxn_vector_matrix, model, targetRxn, rxnListInput, isGeneList)
% The fitness function
%
% USAGE:
%
%    [val] = optGeneFitnessTilt(rxn_vector_matrix, model, targetRxn, rxnListInput, isGeneList)
%
% INPUTS:
%    rxn_vector_matrix:    reactions vectors in a matrix
%    model:                model structure
%    targetRxn:            target reactions
%    rxnListInput:         list of reactions
%    isGeneList:           bolean checking if it is a gene list
%
% OUTPUT:
%    val:                  fitness value

global MaxKnockOuts
global CBT_LP_SOLVER
%size(rxn_vector_matrix)

%rxnGeneMat is a required field for this function, so if it does not exist,
%build it.
if isGeneList && ~isfield(model,'rxnGeneMat')
    model = buildRxnGeneMat(model);
end

popsize = size(rxn_vector_matrix,1);
val = zeros(1,popsize);

for i = 1:popsize
    rxn_vector = rxn_vector_matrix(i,:);
    rxnList = rxnListInput(logical(rxn_vector));

    %see if we've done this before
    val_temp = memoize(rxn_vector);
    if ~ isempty(val_temp)
        val(i) = val_temp;
        continue;
    end

    % check to see if mutations is above the max number allowed
    nummutations = sum(rxn_vector);
    if nummutations > MaxKnockOuts
        continue;
    end

    % generate knockout.
    if isGeneList
        modelKO = deleteModelGenes(model, rxnList);
    else % is reaction list
        [isValidRxn,removeInd] = ismember(rxnList,model.rxns);
        removeInd = removeInd(isValidRxn);
        modelKO = model;
        modelKO.ub(removeInd) = 0;
        modelKO.lb(removeInd) = 0;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % augment BOF (tilt)
    [modelKO] = augmentBOF(modelKO, targetRxn, .001);

    % find growthrate
    if exist('LPBasis', 'var')
        modelKO.LPBasis = LPBasis;
    end

    if ~isempty(CBT_LP_SOLVER) && strcmp(CBT_LP_SOLVER,'cplex')
        [solKO, LPOUT] = solveCobraLPCPLEX(modelKO, 0,1);
        LPBasis = LPOUT.LPBasis;
        growthrate = solKO.obj;
        [~,tar_loc] = ismember(targetRxn,modelKO.rxns);
        minProdAtSetGR = solKO.full(tar_loc);
    else
        solKO = optimizeCbModel(modelKO);
        if isempty(solKO.x)
            continue;
        else
            growthrate = solKO.f;
            [~,tar_loc] = ismember(targetRxn,modelKO.rxns);
            minProdAtSetGR = solKO.x(tar_loc);
        end
    end

    % check to ensure that GR is above a certain value
    if growthrate < .10
        continue;
    end

    % %    display('second optimization');
    %     % find the lowesest possible production rate (a hopefully high number)
    %     % at the max growth rate minus some set factor gamma (a growth rate slightly
    %     % smaller than the max). A positive value will eliminate solutions where the
    %     % production envelope has a vertical line at the max GR, a "non-unique"
    %     % solution. Set value to zero if "non-unique" solutions are not an issue.
    %     gamma = 0.01; % proportional to Grwoth Rate (hr-1), a value around 0.5 max.
    %
    %     %find indicies of important vectors
    %     indBOF = find(modelKO.c);
    %     indTar = findRxnIDs(modelKO, targetRxn);
    %     % generate a model with a fixed max KO growth rate
    %     modelKOsetGR = modelKO;
    %     modelKOsetGR.lb(indBOF) = growthrate - gamma; % this growth rate is required as lb.
    %     modelKOsetGR.c = zeros(size(modelKO.c));
    %     modelKOsetGR.c(indTar) = -1; % minimize for this variable b/c we want to look at the very minimum production.
    %
    %     % find the minimum production rate for the targeted reaction.
    %
    % %     solKOsetGR = optimizeCbModel(modelKOsetGR);
    % %     minProdAtSetGR1 = -solKOsetGR.f;  % This should be a negative value b/c of the minimization setup, so -1 is necessary.
    %
    %     if exist('LPBasis2', 'var')
    %         modelKOsetGR.LPBasis = LPBasis2;
    %     end
    %
    %     [solKOsetGR, LPOUT2] = solveCobraLPCPLEX(modelKOsetGR, 0,1);
    %     LPBasis2 = LPOUT2.LPBasis;
    %     minProdAtSetGR = -solKOsetGR.obj;

    % objective function for optGene algorithm = val (needs to be a negative value, since it is
    % a minimization)
    val(i) = -minProdAtSetGR;
    % penalty for a greater number of mutations

    %val(i) = -minProdAtSetGR * (.98^nummutations);

    % select best substrate-specific productivity
    % val(i) = -minProdAtSetGR * (.98^nummutations) * growthrate;

    % check to prevent very small values from being considerered improvments
    if val(i) > -1e-3
        val(i) = 0;
    end

    memoize(rxn_vector, val(i));
end

return;

%% Memoize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% MEMOIZE %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% internal function to speed things up.

function [value] =  memoize(gene_vector, value)
global HTABLE
hashkey = num2str(gene_vector);
hashkey = strrep(hashkey,' ',''); % cut out white space from string (more space efficient).

if nargin == 1
    value = HTABLE.get(hashkey);
    return;
else
    if HTABLE.size() > 50000
        HTABLE = java.util.Hashtable;  %reset the hashtable if more than 50,000 entries.
    end
    HTABLE.put(hashkey, value);
    value = [];
    return;
end

author: opencobra

repo: cobratoolbox


text:
% Dasslc test problem with classical pendulum
% requires files: pend.m

dae_index = 3;		% differential index of the DAE
g   = 9.8;			% gravity acceleration
L	 = 1.0;			% pendulum cord lenght
t0  = 0.0;        % initial value for independent variable
tf  = 10;         % final value for independent variable
y0  = [1 0 0 0 0]'; % initial state variables (overwritten by pend.dat)
rpar=[g L dae_index]; % optional arguments passed to residual and jacobian functions

index = [0 0 0 0 0	% index 0 formulation (with drift-off effect)
		   1 1 1 1 1	% index 1 formulation
			1 1 2 2 2	% index 2 formulation
         1 1 2 2 3];	% index 3 formulation

tspan=[t0:0.001:tf];
[t,y]=dasslc('pend',tspan,y0,rpar,[],[],index(dae_index+1,:),'pend.dat','jacpend');

plot(t,y);

author: Sable

repo: mcbench-benchmarks


text:
hydro = struct();

hydro = readNEMOH(hydro,'../Cylinder/');
% hydro = readWAMIT(hydro,'../../WAMIT/Cylinder/cyl.out',[]);
% hydro = combineBEM(hydro);   % Compare to WAMIT
hydro = radiationIRF(hydro,5,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,5,[],[],[],[]);
% writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
% CONFIGSTATEMACHINE configures the state machine (type of demo, velocity
%                    of the demo, repeat movements, and so on).

%% --- Initialization ---

% If true, the robot CoM will follow a desired reference trajectory (COORDINATOR DEMO ONLY)
Config.LEFT_RIGHT_MOVEMENTS = false;

% If equal to one, the desired values of the center of mass are smoothed internally
Config.SMOOTH_COM_DES       = true;

% If equal to one, the desired values of the postural tasks are smoothed internally
Config.SMOOTH_JOINT_DES     = true;

% Joint torques saturation [Nm]
Sat.torque                  = 60;

% Joint torques rate of change saturation
Sat.uDotMax                 = 300;

% max unsigned difference between two consecutive (measured) joint positions,
% i.e. delta_qj = abs(qj(k) - qj(k-1))
Sat.maxJointsPositionDelta  = 15*pi/180; % [rad]

%% Regularization parameters
Reg.pinvDamp_baseVel        = 1e-7;
Reg.pinvDamp                = 1;
Reg.pinvTol                 = 1e-5;
Reg.KP_postural             = 0.1;
Reg.KD_postural             = 0;
Reg.HessianQP               = 1e-7;

%% State Machine configuration

% time between two yoga positions
StateMachine.joints_pauseBetweenYogaMoves = 5;

% contact forces threshold
StateMachine.wrench_thresholdContactOn    = 50;
StateMachine.wrench_thresholdContactOff   = 100;

% threshold on CoM and joints error
StateMachine.CoM_threshold                = 0.01;
StateMachine.joints_thresholdNotInContact = 5;
StateMachine.joints_thresholdInContact    = 50;

% initial state for state machine
StateMachine.initialState                 = 1;

% other configuration parameters for state machine
StateMachine.tBalancing                   = 1;
StateMachine.tBalancingBeforeYoga         = 1;
StateMachine.yogaExtended                 = true;
StateMachine.skipYoga                     = false;
StateMachine.demoOnlyBalancing            = false;
StateMachine.demoStartsOnRightSupport     = false; % If false, the Yoga demo is performed on the left foot first
StateMachine.yogaAlsoOnRightFoot          = false; % TO DO: yoga on both feet starting from right foot (not available for now)

%%%% List of possible "Yoga in loop" %%%%

% the robot will repeat the FULL DEMO (two feet balancing, yoga on left
% foot, back on two feet, yoga right foot, back on two feet). The demo is
% repeated until the user stops the Simulink model. This option is ignored
% if Sm.demoStartsOnRightSupport = true.
StateMachine.twoFeetYogaInLoop            = false;

% the robot will repeat the ONE FOOT yoga for the number of times the user
% specifies in the Sm.yogaCounter option. The robot WILL NOT go back to two
% feet balancing in between to consecutive yoga. WARNING: if the option
% Sm.yogaAlsoOnRightFoot is true, then the robot will repeat first the yoga
% on left foot for the number of times the user specifies in the Sm.yogaCounter,
% and then it will repeat the yoga on the right foot for the same number of times.
StateMachine.oneFootYogaInLoop            = false;
StateMachine.yogaCounter                  = 5;

%% Parameters for motors reflected inertia

% transmission ratio (1/N)
Config.Gamma                 = 0.01*eye(ROBOT_DOF);

% modify the value of the transmission ratio for the hip pitch.
% TODO: avoid to hard-code the joint numbering
Config.Gamma(end-5, end-5)   = 0.0067;
Config.Gamma(end-11,end-11)  = 0.0067;

% motors inertia (Kg*m^2)
legsMotors_I_m               = 0.0827*1e-4;
torsoPitchRollMotors_I_m     = 0.0827*1e-4;
torsoYawMotors_I_m           = 0.0585*1e-4;
armsMotors_I_m               = 0.0585*1e-4;

% add harmonic drives reflected inertia
if Config.INCLUDE_HARMONIC_DRIVE_INERTIA

    legsMotors_I_m           = legsMotors_I_m + 0.054*1e-4;
    torsoPitchRollMotors_I_m = torsoPitchRollMotors_I_m + 0.054*1e-4;
    torsoYawMotors_I_m       = torsoYawMotors_I_m + 0.021*1e-4;
    armsMotors_I_m           = armsMotors_I_m + 0.021*1e-4;
end

Config.I_m                   = diag([torsoPitchRollMotors_I_m*ones(2,1);
                                     torsoYawMotors_I_m;
                                     armsMotors_I_m*ones(8,1);
                                     legsMotors_I_m*ones(12,1)]);

% parameters for coupling matrices. Updated according to the wiki:
%
% http://wiki.icub.org/wiki/ICub_coupled_joints
%
% and corrected according to https://github.com/robotology/robots-configuration/issues/39
t            = 0.615;
r            = 0.022;
R            = 0.04;

% coupling matrices
T_LShoulder  = [-1  0  0;
                -1 -t  0;
                 0  t -t];

T_RShoulder  = [ 1  0  0;
                 1  t  0;
                 0 -t  t];

T_torso      = [ 0.5    -0.5     0;
                 0.5     0.5     0;
                 r/(2*R) r/(2*R) r/R];

if Config.INCLUDE_COUPLING

    Config.T = blkdiag(T_torso,T_LShoulder,1,T_RShoulder,1,eye(12));
else
    Config.T = eye(ROBOT_DOF);
end

% gain for feedforward term in joint torques calculation. Valid range: a
% value between 0 and 1
Config.K_ff  = 0;

% Config.USE_DES_JOINT_ACC_FOR_MOTORS_INERTIA if true, the desired joints
% accelerations are used for computing the feedforward term in joint
% torques calculations. Not effective if Config.K_ff = 0.
Config.USE_DES_JOINT_ACC_FOR_MOTORS_INERTIA = false;

%% Constraints for QP for balancing

% The friction cone is approximated by using linear interpolation of the circle.
% So, numberOfPoints defines the number of points used to interpolate the circle
% in each cicle's quadrant
numberOfPoints               = 4;
forceFrictionCoefficient     = 1/3;
torsionalFrictionCoefficient = 1/75;
fZmin                        = 10;

% physical size of the foot
feet_size                    = [-0.07  0.12  ;    % xMin, xMax
                                -0.045 0.045 ];   % yMin, yMax

% Compute contact constraints (friction cone, unilateral constraints)
[ConstraintsMatrix, bVectorConstraints] = wbc.computeRigidContactConstraints ...
    (forceFrictionCoefficient, numberOfPoints, torsionalFrictionCoefficient, feet_size, fZmin);

author: robotology

repo: whole-body-controllers


text:
clear;clc;
ts = 0 : 550;   % 时间

x0=[0.000001];   % 最初感染率，百万分之一

% lambda   日接触率
lambdas = [0.15, 0.07, 0.05, 0.025, 0.0125];
infective_matrix = zeros(length(ts),length(lambdas));

for i = 1 : length(lambdas)
  lambda = lambdas(i);
  [t,x] = ode45(@(t,y) ill_si(t,y,lambda), ts, x0);
  infective_matrix(:,i) = x(:,1);
end

plot(t, infective_matrix);
legend(num2str(lambdas(1)), num2str(lambdas(2)), num2str(lambdas(3)), num2str(lambdas(4)), num2str(lambdas(5)));
title('SI模型');
grid;

repo: NumericalAnalysis


text:
function [saveampname]=ps_load_mean_amp()
%PS_LOAD_MEAN_AMP Load mean_amp file into a matlab workspace
%
%   Andy Hooper, June 2006
%
%   ======================================================================
%   02/2010 AH: Reduce memory needs by reading in chunks
%   11/2010 AH: Do amplitude merge here instead of ps_merge_patches.m
%   11/2010 MA: Fix when length(amp_bit) returns zero at the end
%   ======================================================================

if exist('patch.in','file')
    patch=load('patch.in');
    width=patch(2)-patch(1)+1;
else
    widthname='width.txt';
    if ~exist(widthname,'file')
        widthname= ['../',widthname];
    end
    width=load(widthname);
end

ampname=['./mean_amp.flt'];
saveampname=['./amp_mean.mat'];

if ~exist(ampname,'file')
    ampname=['../mean_amp.flt'];
    saveampname=['../amp_mean.mat'];
end

if ~exist(ampname,'file')
  ampname=['./mean_amp.flt'];
  saveampname=['./amp_mean.mat'];
  fprintf('   Merging mean amplitude files\n')
  widthname='width.txt';
  if ~exist(widthname,'file')
      widthname= ['../',widthname];
  end
  width=load(widthname);
  amp=zeros(width,1,'single');

  if exist('./patch.list','file')
    dirname=struct;
    fid=fopen('patch.list','r');
    i=0;
    while feof(fid)==0
        i=i+1;
        dirname(i).name=fgetl(fid);
    end
    fclose(fid);
  else
    dirname=dir('PATCH_*');
  end
  n_patch=length(dirname);
  for i=1:n_patch
   if ~isempty(dirname(i).name)
    cd(dirname(i).name);
    if exist('./mean_amp.flt','file')
        pxy=load('patch.in');
        fid=fopen('mean_amp.flt');
        amp(pxy(1):pxy(2),pxy(3):pxy(4))=fread(fid,[pxy(2)-pxy(1)+1,inf],'float=>single');
        fclose(fid);
    end
    cd ..
   end
  end

  fid=fopen('./mean_amp.flt','w');
  fwrite(fid,amp,'float');
  fclose(fid);
  clear amp

%error('Cannot find mean_amp.flt in this directory or parent')
end

amp_high=[];
amp_low=[];
fid=fopen(ampname,'r');
while ~feof(fid) % first read through
  amp_bit=fread(fid,[width,1000],'float');
  amp_bit(isnan(amp_bit))=1;
  amp_bit=sort(amp_bit(:));
  if length(amp_bit)==0  % MA fix
    break                % length(amp_bit) returns zero at the end
  end
  amp_high=[amp_high;amp_bit(round(length(amp_bit)*0.91):end)];
  amp_low=[amp_low;amp_bit(round(length(amp_bit)*0.01):end)];
end
fclose(fid);

amp_low=min(amp_low);
amp_high=sort(amp_high);
amp_min=min(amp_low);
amp_max=log(amp_high(ceil(length(amp_high)/2))-amp_min);
clear amp_high amp_low

fid=fopen(ampname,'r');
amp_mean=[];
while ~feof(fid) % second read through
  amp_bit=fread(fid,[width,1000],'float');
  amp_bit=amp_bit'-amp_min;
  amp_bit(amp_bit<1)=1;
  amp_bit=log(amp_bit);
  amp_bit(isnan(amp_bit))=0;
  amp_bit(amp_bit>amp_max)=amp_max;
  amp_bit=uint8(amp_bit*255/amp_max);
  amp_mean=[amp_mean;amp_bit];
  if length(amp_bit)==0  % MA fix
  break
  end
end

save(saveampname,'amp_mean');

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
% RENORM_PARENT_2D Renormalization for 2d scattering
%
% Usage
%   Sx_rn = renorm_parent_2d(Sx)
%
% Input
%   Sx (cell): output of 2d scattering
%
% Output
%   Sx_rn (cell): the renormalized 2d scattering
%
% Description
%   This function will renormalize every 2nd order node by its parent
%   (hence his name renorm_PARENT_2d)
%
% See also
%   SCAT, REMORM_PARENT_3D

function Sx_rn = renorm_parent_2d(Sx)

    Sx_rn = Sx;

    % a function handle to find the parent node
    parent = @(p)(find(Sx{2}.meta.j(1,:) == Sx{3}.meta.j(1,p) &...
        Sx{2}.meta.theta(1,:) == Sx{3}.meta.theta(1,p) & ...
        Sx{2}.meta.q(1,:) == Sx{3}.meta.q(1,p)));

    % for each signal in order 2, divide by its ancestor
    for p = 1:numel(Sx{3}.signal)
        Sx_rn{3}.signal{p} = Sx_rn{3}.signal{p}./Sx{2}.signal{parent(p)};
    end

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet