text:
clc
clear
close all

% To run this script we recommend running example2 first to create all the
% necessary outpouts from GLMsingle that are going to be reused here.

% This script shows how to find single-trial beta weights in the output of
% the GLMsingle. We will show how to average them to create one response
% to each condition rather then a response to each trial. This will produce
% one beta weight for each condition. Additionaly we will show how to
% calculatea t-statistic for each condition and a contrast between two
% example conditions in the fLoc experiment (number vs. face) using single
% trial betas.

% load design
load('./data/nsdflocexampledataset.mat')

figure(1);clf

%Show design matricies.
for d = 1:length(design)
    subplot(2,2,d)
    imagesc(design{d}); colormap gray; drawnow
    xlabel('Conditions')
    ylabel('TRs')
    title(sprintf('Design matrix for run%i',d))
end

% There are 10 differenet conditions in this dataset.

%%
% This NSD fLOC scan session has 6 repetitions of each condition per run.
% In the code below, we are attempting to locate the indices in the beta
% weight GLMsingle outputs modelmd(x,y,z,trials) that correspond to
% repeated conditions.

% Consolidate design matrices
designALL = cat(1,design{:});

% Construct a vector containing 1-indexed condition numbers in
% chronological order.

corder = [];
for p=1:size(designALL,1)
    if any(designALL(p,:))
        corder = [corder find(designALL(p,:))];
    end
end

% We will now find indices for each condition. First lets find all unique
% condition in the corder list

condition_list = unique(corder);

% We will now create a structure where each cell containst indices for
% different condition.

condition_ind = cell(length(condition_list),1);
for i = 1 : length(condition_list)

    condition_ind{i} =  find(corder == condition_list(i));

end

% Knowing the indicies of each condition we can now load an example output
% of GLMsingle and create 1 beta weight for each condition. In the example
% 1 there are 10 unique conditions.

results = load('./example2outputs/GLMsingle/TYPED_FITHRF_GLMDENOISE_RR.mat');
betas = results.modelmd;
sz = size(betas);
betas_average = zeros([sz(1) sz(2) sz(3) length(condition_ind),1]);

for i = 1 : length(condition_list)

    betas_average(:,:,:,i) = nanmean(betas(:,:,:,condition_ind{i}),4);

end

% Initially there were 240 betas (each of 10 conditions was repeated 24
% times.

fprintf('GLMsingle output had %i betas before averaging \n',size(betas,4))

% After averaging there is 1 beta weight for each condition

fprintf('GLMsingle output has %i betas after averaging \n',size(betas_average,4))

%% Calculate t-statistic for each condition
% We calculate t-stat as the mean over the standard deviation across all
% repetitions of the same condition.

t_stat =  zeros([sz(1) sz(2) sz(3) length(condition_ind),1]);

for n = 1 : length(condition_list)

    t_stat(:,:,:,n) = nanmean(betas(:,:,:,condition_ind{i}),4) ./ (std(betas(:,:,:,condition_ind{i}),[],4)./sqrt(length(condition_ind{i})));

end

%% Calculate a contrast between condition 2 (number) and condition 5 (adult face)

[~,p,~,stats] = ttest2(betas(:,:,:,condition_ind{2}),betas(:,:,:,condition_ind{5}),'dim',4);

%% Plot estiamted contrasts
slices = [20 10];

underlay = data{1}(:,:,:,1);

things2plot = {'t_stat';'stats.tstat'};

figure(1);clf
for f = 1 : length(things2plot)

    slice = slices(f);
    subplot(1,2,f)
    overlay  = eval(things2plot{f});
    overlay  = overlay(:,:,slice);

    brainmask = (overlay < -3 | overlay > 3);

    underlay_im = cmaplookup(underlay,min(underlay(:)),max(underlay(:)),[],gray(256));
    overlay_im  = cmaplookup(overlay,-4,4,[],cmapsign2);

    hold on
    imagesc(squeeze(underlay_im(:,:,slice,:)));
    imagesc(overlay_im,'AlphaData',brainmask);
    axis image
    axis off
    set(gca,'FontSize',14)
    title(sprintf('%s, slice = %i',things2plot{f},slice),'Interpreter','None')

    colormap(cmapsign2)
    c = colorbar;
    c.Ticks = [0 0.5 1];
    c.TickLabels = {'-3';'0';'3'};

end

% The t-statistic (t_stat) shows high values in the visual cortex while the
% t-test between two conditions shows positive voxels in more ventral and
% lateral portions of the brain.

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function y = BF_PreProcess(y,preProcessHow)
% BF_PreProcess   Preprocess a time series, y

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

if ~isempty(preProcessHow)
    switch preProcessHow
    case 'diff1'
        % Takes incremental differences of the input time series
        y = diff(y);
    case 'rescale_tau'
        % Coarse-graining at a given scale, as in multiscale entropy approaches
        % Find first zero of the autocorrelation function
        tau = CO_FirstCrossing(y,'ac',0,'discrete');
        % Buffer the time series into nonoverlapping windows of length tau
        y_buffer = BF_MakeBuffer(y,tau);
        % Mean each window to get a coarse-grained time series
        y = mean(y_buffer,2);
    otherwise
        error('Unknown preprocessing setting: ''%s''',preProcessHow);
    end
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [version released] = manopt_version()
% Returns the version of the Manopt package you are running, as a vector.
%
% function [version released] = manopt_version()
%
% version(1) is the primary version number.
% released is the date this version was released, in the same format as the
% date() function in Matlab.

    version = [1, 0, 6];
    released = '25-Jun-2014';

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function make_noise_model
%% Input params about your paths and data
%start by building a rez structure -- will be used to call the whitening
%and filtering functions in KS2

%these need to match the current location of KS2 and npy

addpath(genpath('D:\KS2\Kilosort2-master')) % path to kilosort folder
addpath('D:\KS2\npy-matlab-master\')

% params normally read in from a config file. These need to match the model
% data set from which the noise will be extracted.

% sample rate
ops.fs = 30000;
ops.bitPerUV = 0.42667;

% chanMap for your input data set. Make sure it includes reference channels
% where they occur in the data
ops.chanMap = 'D:\SC026_080819_g0_Imec2\SC026_noise_model\neuropixPhase3B1_kilosortChanMap.mat';
[~,chanMapName,~] = fileparts(ops.chanMap);

%path to your data file; results will also be written to this folder
rootZ = 'D:\SC026_080819_g0_Imec2\SC026_noise_model\';
%name of the binary file
dataName = 'SC026_080819_g0_tcat.imec2.ap.bin';

%directory for temporary whitened data filt
rootH = 'D:\KS2\kilosort_datatemp\';
ops.trange = [0 inf]; % time range to use when extracting the whitening matrix
ops.NchanTOT    = 385; % total number of channels in your recording, including digital
%Time to sample in sec, need to sample at least enough to freq down to the high pass limit
clipDur = 0.25;
%%

% frequency for high pass filtering (150)
ops.fshigh = 150;

%processing params usually read in from config. These should not be changed
ops.GPU                 = 1; % has to be 1, no CPU version yet, sorry
ops.useRAM              = 0;
ops.ntbuff              = 64;    % samples of symmetrical buffer for whitening and spike detection
ops.NT                  = 64*1024+ ops.ntbuff; % must be multiple of 32 + ntbuff. This is the batch size (try decreasing if out of memory).
ops.whiteningRange      = 32; % number of channels to use for whitening each channel
ops.nSkipCov            = 25; % compute whitening matrix from every N-th batch
ops.scaleproc           = 100;   % int16 scaling of whitened data

ops.fproc       = fullfile(rootH, 'temp_wh.imec.ap.bin'); % proc file on a fast SSD

ops.fbinary = fullfile(rootZ, dataName);
ops.minfr_goodchannels = 0; %don't remove any channels due to low firing rate

% call preprocessDataSub to filter data and create the whitening matrix.
ops.trange = [0 60]; % time range to use when extracting the whitening matrix
[rez, ops] = preprocessDataForNoise(ops);

%calculate fft for each channel in the whitened data (stored in ops.fproc)
tic;
[fftSamp, nT_fft] = sampFFT( rez, clipDur );
fprintf( 'time to calc fft: %f\n', toc );
%use the fft to generate noise for each channel for a time period
%Can only make batches of nT points at a time

    noiseT = 2.2;           %in seconds, chosen for ~66000 points
    noiseSamp = noiseT*rez.ops.fs;
    nChan = numel(rez.ops.chanMap);
    eNoise = zeros( noiseSamp, nChan, 'single' );

    noiseBatch = ceil(noiseSamp/nT_fft);
    lastWind = noiseSamp - (noiseBatch-1)*nT_fft; %in samples
    fprintf( 'noiseSamp, batchSamp, lastWind: %d, %d, %d\n', noiseSamp, nT_fft, lastWind);

    for j = 1:nChan
            for i = 1:noiseBatch-1
                tStart = (i-1)*nT_fft+1;
                tEnd = i * nT_fft;
                eNoise(tStart:tEnd,j) = fftnoise(fftSamp(:,j),1);
            end
            %for last batch, call one more time and truncate
            tempNoise = fftnoise(fftSamp(:,j),1);
            tStart = (noiseBatch-1)*nT_fft+1;
            tEnd = noiseSamp;
            eNoise(tStart:tEnd,j) = tempNoise(1:lastWind);
    end

    selectChan = [10,11,12,13,14];

    tR = [4200:5200]; %just looking at 1000 samples
    h = figure(1);
    for k = 1:numel(selectChan)
        currNoise = eNoise(tR,selectChan(k));
        plot( tR, currNoise + k*500 );
        hold on;
    end
    title("1000 samples of generated noise, before unwhitening");
    hold off

    %unwhiten this noise array
    eNoise_unwh = eNoise/rez.Wrot;

    %to get the final noise array, map to an array including all channels
    nAllChan = max(rez.ops.chanMap);
    eNoise_final = zeros(noiseSamp, nAllChan, 'single');
    eNoise_final(:,rez.ops.chanMap) = eNoise_unwh;

    tR = [4200:5200]; %plot 1000 samples post unwhitening
    h = figure(2);
    for k = 1:numel(selectChan)
        currNoise = eNoise_final(tR,selectChan(k));
        plot( tR, currNoise + k*50 );
        hold on;
    end
    title('1000 samples of generated noise, after unwhitening');
    hold off

    nm.chanMapName = chanMapName;
    nm.fft = fftSamp;
    nm.nt = nT_fft;
    nm.Wrot = rez.Wrot;
    nm.chanMap = rez.ops.chanMap;
    nm.bitPerUV =  rez.ops.bitPerUV;

%save rez file
fprintf('Saving rez file  \n');
fname = fullfile(rootZ, 'rez.mat');
save(fname, 'rez', '-v7.3');

%save noise model file
fprintf('Saving nm file  \n');
fname = fullfile(rootZ, 'noiseModel.mat');
save(fname, 'nm', '-v7.3');

end

function [rez, ops] = preprocessDataForNoise(ops)

%build rez structure -- code taken from preprocessDataSub
%Only changes are to skip checking for low spike rate samples and
%write out the data as rows = time -- just a convenience so it
%it can be read in using the same code

ops.nt0 	= getOr(ops, {'nt0'}, 61);
ops.nt0min  = getOr(ops, 'nt0min', ceil(20 * ops.nt0/61));

NT       = ops.NT ;
NchanTOT = ops.NchanTOT;

bytes = get_file_size(ops.fbinary);
nTimepoints = floor(bytes/NchanTOT/2);
ops.tstart = ceil(ops.trange(1) * ops.fs);
ops.tend   = min(nTimepoints, ceil(ops.trange(2) * ops.fs));
ops.sampsToRead = ops.tend-ops.tstart;
ops.twind = ops.tstart * NchanTOT*2;

Nbatch      = ceil(ops.sampsToRead /(NT-ops.ntbuff));
ops.Nbatch = Nbatch;

[chanMap, xc, yc, kcoords, NchanTOTdefault] = loadChanMap(ops.chanMap);
ops.NchanTOT = getOr(ops, 'NchanTOT', NchanTOTdefault);
ops.igood = true(size(chanMap));

ops.Nchan = numel(chanMap);
ops.Nfilt = getOr(ops, 'nfilt_factor', 4) * ops.Nchan;

rez.ops         = ops;
rez.xc = xc;
rez.yc = yc;

rez.xcoords = xc;
rez.ycoords = yc;

rez.ops.chanMap = chanMap;
rez.ops.kcoords = kcoords;

NTbuff      = NT + 4*ops.ntbuff;

% by how many bytes to offset all the batches
rez.ops.Nbatch = Nbatch;
rez.ops.NTbuff = NTbuff;
rez.ops.chanMap = chanMap;

%build whitening matrix. This function filters before calculating the
%cross correlation

tic;

fprintf('Time %3.0fs. Computing whitening matrix.. \n', toc);

Wrot = get_whitening_matrix(rez);

%Apply the whitening matrix to a subset of the original data

ops.tstart = ceil(ops.trange(1) * ops.fs);
ops.tend   = min(nTimepoints, ceil(ops.trange(2) * ops.fs));
ops.sampsToRead = ops.tend-ops.tstart;
ops.twind = ops.tstart * NchanTOT*2;
Nbatch      = ceil(ops.sampsToRead /(NT-ops.ntbuff));
ops.Nbatch = Nbatch;

fprintf('Time %3.0fs. Loading raw data and applying filters... \n', toc);

fid         = fopen(ops.fbinary, 'r');
if ~ops.useRAM
    fidW        = fopen(ops.fproc,   'w');
    DATA = [];
else
    DATA = zeros(NT, rez.ops.Nchan, Nbatch, 'int16');
end
% load data into batches, filter, compute covariance
if isfield(ops,'fslow')&&ops.fslow<ops.fs/2
    [b1, a1] = butter(3, [ops.fshigh/ops.fs,ops.fslow/ops.fs]*2, 'bandpass');
elseif isfield(ops,'fshigh')
    [b1, a1] = butter(3, ops.fshigh/ops.fs*2, 'high');
end

for ibatch = 1:Nbatch
    offset = max(0, ops.twind + 2*NchanTOT*((NT - ops.ntbuff) * (ibatch-1) - 2*ops.ntbuff));
    if offset==0
        ioffset = 0;
    else
        ioffset = ops.ntbuff;
    end
    fseek(fid, offset, 'bof');

    buff = fread(fid, [NchanTOT NTbuff], '*int16');
    if isempty(buff)
        break;
    end
    nsampcurr = size(buff,2);
    if nsampcurr<NTbuff
        buff(:, nsampcurr+1:NTbuff) = repmat(buff(:,nsampcurr), 1, NTbuff-nsampcurr);
    end

    if ops.GPU
        dataRAW = gpuArray(buff);
    else
        dataRAW = buff;
    end
    dataRAW = dataRAW';
    dataRAW = single(dataRAW);
    dataRAW = dataRAW(:, chanMap);

    % subtract the mean from each channel
    dataRAW = dataRAW - mean(dataRAW, 1);

    if (isfield(ops,'fshigh'))
        datr = filter(b1, a1, dataRAW);
        datr = flipud(datr);
        datr = filter(b1, a1, datr);
        datr = flipud(datr);
    else
        datr = dataRAW;
    end

    % CAR, common average referencing by median
    if getOr(ops, 'CAR', 1)
        datr = datr - median(datr, 2);
    end

    datr = datr(ioffset + (1:NT),:);

    %apply whitening
    datr    = datr * Wrot;

    if ops.useRAM
        DATA(:,:,ibatch) = gather_try(datr);
    else
        datcpu  = gather_try(int16(datr));
        %next step will be reading in this file to get FTs
        %Most convenient for that process to have the data in the original
        %format (row = time, column = channel index)
        fwrite(fidW, datcpu', 'int16');
    end
end

Wrot        = gather_try(Wrot);
rez.Wrot    = Wrot;

fclose(fidW);
fclose(fid);

end

function [fftSamp, nT] = sampFFT( rez, clipDur )

    ops = rez.ops;

    %get file of whitened data
    datFile = ops.fproc;

    %read in some data
    nChansInFile = numel(ops.chanMap);  % channels in whitenened data, excludes ref chans in original
    d = dir(ops.fproc);
    nSamps = d.bytes/2/nChansInFile;
    tSkip = 1.0; % skip first tSkip seconds

    sampStart = round(ops.fs*tSkip);
    nClipSamps = round(ops.fs*clipDur);
    mmf = memmapfile(datFile, 'Format', {'int16', [nChansInFile nSamps], 'x'});
    thisDat = (double(mmf.Data.x(:, (1:nClipSamps)+sampStart)));
    %subtract the DC
    thisDat = bsxfun(@minus, thisDat, mean(thisDat,2));
    [~,nT] = size(thisDat);
    %data is formatted as Nchan rows by nt columns. fft returns the
    %transform of the columns, so transpose
    fftSamp = fft(thisDat');

end

function noise=fftnoise(f,Nseries)
% Generate noise with a given power spectrum.
% Useful helper function for Monte Carlo null-hypothesis tests and confidence interval estimation.
%
% noise=fftnoise(f[,Nseries])
%
% INPUTS:
% f: the fft of a time series (must be a column vector)
% Nseries: number of noise series to generate. (default=1)
%
% OUTPUT:
% noise: surrogate series with same power spectrum as f. (each column is a surrogate).
%
%   --- Aslak Grinsted (2009)
if nargin<2
    Nseries=1;
end
f=f(:);     %ensures f is a column vector
N=length(f);
Np=floor((N-1)/2);
phases=rand(Np,Nseries)*2*pi;
phases=complex(cos(phases),sin(phases)); % this was the fastest alternative in my tests.
f=repmat(f,1,Nseries);
f(2:Np+1,:)=f(2:Np+1,:).*phases;
f(end:-1:end-Np+1,:)=conj(f(2:Np+1,:));
noise=real(ifft(f,[],1));

end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
% UNPAD_SIGNAL Remove de padding from PAD_SIGNAL
%
% Usage
%    x = UNPAD_SIGNAL(y, resolution, target_sz, center)
%
% Input
%    y (numeric): The signal to be unpadded.
%    resolution (int): The resolution of the signal (as a power of 2), with
%        respect to the original, unpadded version.
%    target_sz (numeric): The size of the original, unpadded version. Combined
%        with resolution, the size of the output y is given by
%        target_sz.*2.*(-resolution).
%    center (boolean, optional): If true, extracts the center part of y, oth-
%        erwise extracts the (upper) left corner (default false).
%
% Output
%    x (numeric): The extracted unpadded signal
%
% Description
%    To handle boundary conditions, a signal is often padded using PAD_SIGNAL
%    before being convolved with CONV_SUB_1D or CONV_SUB_2D. After this, the
%    padding needs to be removed to recover a regular signal. This is achieved
%    using UNPAD_SIGNAL, which takes the padded, convolved signal y as input,
%    as well as its resolution relative to the original, unpadded version,
%    and the size of this original version. Using this, it extracts the
%    coefficients in y that correspond to the domain of the original signal.
%    If the center flag was specified during PAD_SIGNAL, it is specified here
%    again in order to extract the correct part.
%
% See Also
%    PAD_SIGNAL
function x = unpad_signal(x, res, target_sz, center)
    if nargin < 4
        center = 0;
    end

    padded_sz = size(x);

    padded_sz = padded_sz(1:length(target_sz));

    offset = 0.*target_sz;

    if center
        offset = (padded_sz.*2.^res-target_sz)/2;
    end

    offset_ds = floor(offset./2.^res);
    target_sz_ds = 1+floor((target_sz-1)./2.^res);

    switch length(target_sz)
        case 1
            x = x(offset_ds + (1:target_sz_ds),:,:);
        case 2
            x = x(offset_ds(1) + (1:target_sz_ds(1)), ...
                offset_ds(2) + (1:target_sz_ds(2)), :);
    end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [LLF, likelihoods, errors] = sarimax_likelihood(parameters,p,q,constant,seasonal,y,x,sigma)
% PURPOSE:
%   Likelihood function for armaxfilter
%
% USAGE:
%   [LLF, LIKELIHOODS, ERRORS] = armaxfilter_likelihood(PARAMETERS,P,Q,CONSTANT,Y,X,M)
%   [LLF, LIKELIHOODS, ERRORS] = armaxfilter_likelihood(PARAMETERS,P,Q,CONSTANT,Y,X,M)
%
% INPUTS:
%   PARAMETERS  - A vector of GARCH process aprams of the form [constant, arch, garch]
%   P           - Vector containing lag indices of the AR component
%   Q           - Vector containing lag indices of the MA component
%   CONSTANT    - Value indicating whether the model contains a constant (1) or not (0)
%   Y           -
%   X           -
%   M           - Index to first element to use in the recursive residual calculation
%   SIGMA       - T by 1 vector
%
% OUTPUTS:
%   LLF         - Minus 1 times the log likelihood
%   LIKELIHOODS - Time series of likelihoods
%   ERRORS      - Time series of model errors
%
% COMMENTS:
%
% See also armaerrors

% Author: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 4/1/2004

m = size(x,2);
[armaParameters, p, q] = sarma2arma(parameters(m+1:length(parameters)), p, q, seasonal);
parameters = [parameters(1:m);armaParameters];

e = armaxerrors(parameters,p,q,constant,y,x,m,ones(size(y)));
T = length(e);

% Do not divide e by sigma since this is done in armaxerrors
likelihoods =  0.5*(2*log(sigma) + (e./sigma).^2 + log(2*pi));

if isempty(p)
    p = 0;
end
if isempty(q)
    q = 0;
end
if max(q)>max(p) % prepend y and x, if needed
    t = (max(q)-max(p))+1:T;
else
    t = 1:T;
end

likelihoods = likelihoods(t);
errors=e(t);
LLF = sum(likelihoods);

if isnan(LLF)
    LLF=1e7;
end

author: bashtage

repo: mfe-toolbox


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

    num_images = length(image_roidb);
    if conf.bbox_class_agnostic
        num_classes = 1;
    else
        % Infer number of classes from the number of columns in gt_overlaps
        num_classes = size(image_roidb(1).overlap, 2);
    end
    % Sample random scales to use for each image in this batch
    random_scale_inds = randi(length(conf.scales), num_images, 1);

    if conf.batch_size > 0
        assert(mod(conf.batch_size, num_images) == 0, ...
            sprintf('num_images %d must divide BATCH_SIZE %d', num_images, conf.batch_size));

        rois_per_image = conf.batch_size / num_images;
        fg_rois_per_image = round(rois_per_image * conf.fg_fraction);
    else
        rois_per_image = inf;
        fg_rois_per_image = inf;
    end

    % Get the input image blob
    [im_blob, im_scales] = get_image_blob(conf, image_roidb, random_scale_inds);

    % build the region of interest and label blobs
    rois_blob = zeros(0, 5, 'single');
    labels_blob = zeros(0, 1, 'single');
    bbox_targets_blob = zeros(0, 4 * (num_classes+1), 'single');
    bbox_loss_blob = zeros(size(bbox_targets_blob), 'single');

    for i = 1:num_images
        [labels, ~, im_rois, bbox_targets, bbox_loss] = ...
            sample_rois(conf, image_roidb(i), fg_rois_per_image, rois_per_image);

        % Add to ROIs blob
        feat_rois = rfcn_map_im_rois_to_feat_rois(conf, im_rois, im_scales(i));
        batch_ind = i * ones(size(feat_rois, 1), 1);
        rois_blob_this_image = [batch_ind, feat_rois];
        rois_blob = [rois_blob; rois_blob_this_image];

        % Add to labels, bbox targets, and bbox loss blobs
        labels_blob = [labels_blob; labels];
        bbox_targets_blob = [bbox_targets_blob; bbox_targets];
        bbox_loss_blob = [bbox_loss_blob; bbox_loss];
    end

    % permute data into caffe c++ memory, thus [num, channels, height, width]
    im_blob = im_blob(:, :, [3, 2, 1], :); % from rgb to brg
    im_blob = single(permute(im_blob, [2, 1, 3, 4]));
    rois_blob = rois_blob - 1; % to c's index (start from 0)
    rois_blob = single(permute(rois_blob, [3, 4, 2, 1]));
    labels_blob = single(permute(labels_blob, [3, 4, 2, 1]));
    bbox_targets_blob = single(permute(bbox_targets_blob, [3, 4, 2, 1]));
    bbox_loss_blob = single(permute(bbox_loss_blob, [3, 4, 2, 1]));

    assert(~isempty(im_blob));
    assert(~isempty(rois_blob));
    assert(~isempty(labels_blob));
    assert(~isempty(bbox_targets_blob));
    assert(~isempty(bbox_loss_blob));

    net_inputs = {im_blob, rois_blob, labels_blob, bbox_targets_blob, bbox_loss_blob};
end

%% Build an input blob from the images in the roidb at the specified scales.
function [im_blob, im_scales] = get_image_blob(conf, images, random_scale_inds)

    num_images = length(images);
    processed_ims = cell(num_images, 1);
    im_scales = nan(num_images, 1);
    for i = 1:num_images
        im = imread(images(i).image_path);
        target_size = conf.scales(random_scale_inds(i));

        [im, im_scale] = prep_im_for_blob(im, conf.image_means, target_size, conf.max_size);

        im_scales(i) = im_scale;
        processed_ims{i} = im;
    end

    im_blob = im_list_to_blob(processed_ims);
end

%% Generate a random sample of ROIs comprising foreground and background examples.
function [labels, overlaps, rois, bbox_targets, bbox_loss_weights] = sample_rois(conf, image_roidb, fg_rois_per_image, rois_per_image)

    [overlaps, labels] = max(image_roidb(1).overlap, [], 2);
%     labels = image_roidb(1).max_classes;
%     overlaps = image_roidb(1).max_overlaps;
    rois = image_roidb(1).boxes;

    % Select foreground ROIs as those with >= FG_THRESH overlap
    fg_inds = find(overlaps >= conf.fg_thresh);
    % Guard against the case when an image has fewer than fg_rois_per_image
    % foreground ROIs
    fg_rois_per_this_image = min(fg_rois_per_image, length(fg_inds));
    % Sample foreground regions without replacement
    if ~isempty(fg_inds)
       fg_inds = fg_inds(randperm(length(fg_inds), fg_rois_per_this_image));
    end

    % Select background ROIs as those within [BG_THRESH_LO, BG_THRESH_HI)
    bg_inds = find(overlaps < conf.bg_thresh_hi & overlaps >= conf.bg_thresh_lo);
    % Compute number of background ROIs to take from this image (guarding
    % against there being fewer than desired)
    bg_rois_per_this_image = rois_per_image - fg_rois_per_this_image;
    bg_rois_per_this_image = min(bg_rois_per_this_image, length(bg_inds));
    % Sample foreground regions without replacement
    if ~isempty(bg_inds)
       bg_inds = bg_inds(randperm(length(bg_inds), bg_rois_per_this_image));
    end
    % The indices that we're selecting (both fg and bg)
    keep_inds = [fg_inds; bg_inds];
    % Select sampled values from various arrays
    labels = labels(keep_inds);
    % Clamp labels for the background ROIs to 0
    labels((fg_rois_per_this_image+1):end) = 0;
    overlaps = overlaps(keep_inds);
    rois = rois(keep_inds, :);

    if conf.bbox_class_agnostic
        assert(all((labels>0) == image_roidb.bbox_targets(keep_inds, 1)));
    else
        assert(all(labels == image_roidb.bbox_targets(keep_inds, 1)));
    end

    % Infer number of classes from the number of columns in gt_overlaps
    num_classes = size(image_roidb(1).overlap, 2);

    [bbox_targets, bbox_loss_weights] = get_bbox_regression_labels(conf, ...
        image_roidb.bbox_targets(keep_inds, :), num_classes);

end

function [bbox_targets, bbox_loss_weights] = get_bbox_regression_labels(conf, bbox_target_data, num_classes)
%% Bounding-box regression targets are stored in a compact form in the roidb.
 % This function expands those targets into the 4-of-4*(num_classes+1) representation used
 % by the network (i.e. only one class has non-zero targets).
 % The loss weights are similarly expanded.
% Return (N, (num_classes+1) * 4, 1, 1) blob of regression targets
% Return (N, (num_classes+1 * 4, 1, 1) blob of loss weights
    if conf.bbox_class_agnostic
        num_classes = 1;
    end

    clss = bbox_target_data(:, 1);
    bbox_targets = zeros(length(clss), 4 * (num_classes+1), 'single');
    bbox_loss_weights = zeros(size(bbox_targets), 'single');
    inds = find(clss > 0);
    for i = 1:length(inds)
       ind = inds(i);
       cls = clss(ind);
       targets_inds = (1+cls*4):((cls+1)*4);
       bbox_targets(ind, targets_inds) = bbox_target_data(ind, 2:end);
       bbox_loss_weights(ind, targets_inds) = 1;
    end
end

repo: R-FCN


text:
% Implement the back propagation of an LSTM layer.
% Author: Xiong Xiao, Temasek Labs, NTU, Singapore.
% Last modified: 13 Oct 2015
%
function [grad, grad_W, grad_b] = B_LSTM(input, LSTM_layer, future_layers)
W = LSTM_layer.W;
if strcmpi(class(input), 'gpuArray'); useGPU=1; else useGPU = 0; end
precision = class(gather(input(1)));

[dim, nFr, nSeg] = size(input);
input = permute(input, [1 3 2]);

nCell = LSTM_layer.dim(1);  % number of LSTM cells in the layer

useHidden = 1;
usePastState = 1;
usePastStateAsFeature = LSTM_layer.usePastState;

future_layer_grad = 0;
for fli = 1:length(future_layers)
    future_layer_grad = future_layer_grad + future_layers{fli}.grad;      % this is the gradient of LSTM output from the cost function
end
grad_ht_cost = future_layer_grad;
grad_ht_cost = permute(grad_ht_cost,[1 3 2]);

ft = LSTM_layer.ft;
it = LSTM_layer.it;
ot = LSTM_layer.ot;
Ct_raw = LSTM_layer.Ct_raw;
Ct = LSTM_layer.Ct;
ht = LSTM_layer.a;
ht = permute(ht, [1 3 2]);
Ct0 = LSTM_layer.Ct0;
ht0 = LSTM_layer.ht0;

% allocate memory for the gradients of gates and states
if useGPU == 0
    grad_xt = zeros(dim, nSeg, nFr, precision);
    grad_ft = zeros(nCell, nSeg, nFr, precision);     % forget gates
    grad_it = zeros(nCell, nSeg, nFr, precision);     % input gates
    grad_ot = zeros(nCell, nSeg, nFr, precision);     % output gates
    grad_Ct_raw = zeros(nCell, nSeg, nFr, precision); % candidate cell states
    grad_Ct = zeros(nCell, nSeg, nFr, precision);     % cell states
    grad_Ct_future = zeros(nCell, nSeg, nFr, precision);     % cell states
    grad_ht = zeros(nCell, nSeg, nFr, precision);     % hidde layer output, i.e. the output of the LSTM layer
    grad_ht_future = zeros(nCell, nSeg, nFr, precision);     % hidde layer output, i.e. the output of the LSTM layer
    grad_zt = zeros(nCell*4,nSeg, nFr, precision);
else
    grad_xt = gpuArray.zeros(dim, nSeg, nFr, precision);
    grad_ft = gpuArray.zeros(nCell, nSeg, nFr, precision);     % forget gates
    grad_it = gpuArray.zeros(nCell, nSeg, nFr, precision);     % input gates
    grad_ot = gpuArray.zeros(nCell, nSeg, nFr, precision);     % output gates
    grad_Ct_raw = gpuArray.zeros(nCell, nSeg, nFr, precision); % candidate cell states
    grad_Ct = gpuArray.zeros(nCell, nSeg, nFr, precision);     % cell states
    grad_Ct_future = gpuArray.zeros(nCell, nSeg, nFr, precision);     % cell states
    grad_ht = gpuArray.zeros(nCell, nSeg, nFr, precision);     % hidde layer output, i.e. the output of the LSTM layer
    grad_ht_future = gpuArray.zeros(nCell, nSeg, nFr, precision);     % hidde layer output, i.e. the output of the LSTM layer
    grad_zt = gpuArray.zeros(nCell*4,nSeg, nFr, precision);
end

for t = nFr:-1:1
    % compute the gradient of ht and Ct that requires gradients from
    % future. At frame nFr, the future gradients are initialized to 0.
    Ct_raw_curr = Ct_raw(:,:,t);
    if useHidden
        grad_ht_curr = grad_ht_future(:,:,t) + grad_ht_cost(:,:,t);
    else
        grad_ht_curr = grad_ht_cost(:,:,t);
    end
    tanh_Ct = tanh(Ct(:,:,t));
    grad_Ct(:,:,t) = grad_Ct_future(:,:,t) + grad_ht_curr .* (1-tanh_Ct.*tanh_Ct) .* ot(:,:,t);

    % compute the gradient of gates and candidate states
    if usePastState
        if t==1
            grad_ft(:,:,t) = grad_Ct(:,:,t) .* Ct0;
            grad_Ct_future0 = grad_Ct(:,:,t) .* ft(:,:,t);
        else
            grad_ft(:,:,t) = grad_Ct(:,:,t) .* Ct(:,:,t-1);
            grad_Ct_future(:,:,t-1) = grad_Ct(:,:,t) .* ft(:,:,t);
        end
    end
    grad_Ct_raw(:,:,t) = grad_Ct(:,:,t) .* it(:,:,t);
    grad_it(:,:,t) = grad_Ct(:,:,t) .* Ct_raw_curr;
    grad_ot(:,:,t) = grad_ht_curr .* tanh(Ct(:,:,t));

    % compute the gradient of the gates before the activation function.

    grad_zCt_raw = grad_Ct_raw(:,:,t) .* (1-Ct_raw_curr.*Ct_raw_curr);
    if 1
        gates = [ft(:,:,t); it(:,:,t); ot(:,:,t)];
        grad_gates = [grad_ft(:,:,t); grad_it(:,:,t); grad_ot(:,:,t)];
        grad_zgates = grad_gates .* gates .* (1-gates);
        grad_zt_curr = [grad_zgates(1:nCell,:); grad_zCt_raw; grad_zgates(nCell+1:end,:)];
    else
        grad_zft = grad_ft(:,:,t) .* ft(:,:,t) .* (1-ft(:,:,t));
        grad_zit = grad_it(:,:,t) .* it(:,:,t) .* (1-it(:,:,t));
        grad_zot = grad_ot(:,:,t) .* ot(:,:,t) .* (1-ot(:,:,t));
        grad_zt_curr = [grad_zft; grad_zCt_raw; grad_zit; grad_zot];
    end
    grad_zt(:,:,t) = grad_zt_curr;

    % compute the gradient of the W, b, and past hidden, past state, and x
    grad_yt = W' * grad_zt_curr;
    if usePastStateAsFeature
        grad_xt(:,:,t) = grad_yt(nCell*2+1:end,:);
    else
        grad_xt(:,:,t) = grad_yt(nCell+1:end,:);
    end
    if t==1
%         grad_ht_future0 = grad_yt(nCell+1:nCell*2);
%         grad_Ct_future0 = grad_Ct_future0 + grad_yt(1:nCell);
%         grad_W = grad_W + grad_zt * [Ct0*usePastStateAsFeature; ht0*useHidden; input(:,t)]';
    else
        grad_ht_future(:,:,t-1) = grad_yt(nCell+1:nCell*2,:);
        if usePastStateAsFeature
            grad_Ct_future(:,:,t-1) = grad_Ct_future(:,:,t-1) + grad_yt(1:nCell,:);
        end
%         grad_W = grad_W + grad_zt * [Ct(:,t-1)*usePastStateAsFeature; ht(:,t-1)*useHidden; input(:,t)]';
    end
%     grad_b = grad_b + grad_zt;
    grad_ht(:,:,t) = grad_ht_curr;
end
if usePastStateAsFeature
    tmpMat = [Ct(:,:,1:nFr-1); ht(:,:,1:nFr-1)*useHidden; input(:,:,2:nFr)];
else
    tmpMat = [ht(:,:,1:nFr-1)*useHidden; input(:,:,2:nFr)];
end
tmpMat = reshape(tmpMat, size(tmpMat,1), nSeg*(nFr-1));
tmpMat2 = reshape(grad_zt(:,:,2:nFr), nCell*4, nSeg*(nFr-1));
grad_W = tmpMat2 * tmpMat';
if usePastStateAsFeature
    grad_W = grad_W + grad_zt(:,:,1) * [Ct0; ht0*useHidden; input(:,:,1)]';
else
    grad_W = grad_W + grad_zt(:,:,1) * [ht0*useHidden; input(:,:,1)]';
end

grad_b = sum(sum(grad_zt,3),2);
grad = permute(grad_xt,[1 3 2]);
end

author: singaxiong

repo: SignalGraph


text:
function addcopyright(organization,folder,copyright_year,update,copyright_string)
%ADDCOPYRIGHT Add copyright info to M-file.
%   ADDCOPYRIGHT(ORGANIZATION) adds the copyright info to all M-files in
%   the current folder as well as subfolders.  The copyright year is the
%   current year as a default.  The organization is set to ORGANIZATION.
%
%   The copyright line that is built up looks like the following:
%
%         %   Copyright <YEAR> <ORGANIZATION>
%
%   For example:
%
%         %   Copyright 2003 The MathWorks, Inc.
%
%   The copyright info is place directly above the Revision line.
%   Specifically is it looking for the following substring:
%
%          %   $Revision
%
%   Where the comment character is placed in the first column.  If it
%   can't find the $Revision keyword (Typically this is right after the
%   help section of the M-file.  See this file as an example), it doesn't
%   know where to insert it.  When it is finished, it will look like the
%   following example:
%
%          %   Copyright 2003 The MathWorks, Inc.
%          %   $Revision...
%
%   ADDCOPYRIGHT(ORGANIZTIONS,FOLDER) adds the copyright info to all
%   M-files in the folder FOLDER and its subfolders.  The copyright year
%   is the current year as a default.  The organization is set to
%   ORGANIZATION.
%
%   ADDCOPYRIGHT(ORGANIZATION,FOLDER,COPYRIGHT_YEAR) adds the copyright
%   info to all M-files in the folder FOLDER and its subfolders.  If
%   FOLDER is [], the default value, the current folder (pwd), is used.
%   The copyright year is set to COPYRIGHT_YEAR.  COPYRIGHT_YEAR may
%   either be a character string (for example, '2003' or '2001-2003') or
%   a numeric.  The organization is set to ORGANIZATION.
%
%   ADDCOPYRIGHT(...,UPDATE) updates the copyright info if UPDATE is set
%   to true (false by default).  If FOLDER is [], the default value, the
%   current folder, is used.  If COPYRIGHT_YEAR is [], the default value,
%   the current year, is used.
%
%   When replacing the copyright info line, it looks for the following string:
%
%          %   Copyright
%
%   Where the comment character is placed in the first column.
%
%   ADDCOPYRIGHT([],FOLDER,[],UPDATE,COPYRIGHT_STRING) ignore the
%   organization and the copyright year to build up a copyright string
%   and use COPYRIGHT_STRING instead.
%
%   Examples:
%   ---------
%   % Supply just the organization name, starting from the current
%   % folder, and use the current year but don't overwrite it if the
%   % copyright info is already there.
%   addcopyright('My Organization')
%
%   % Supply the organization name starting from the "src" folder (found
%   % in the current folder), and use the current year but don't
%   % overwrite it if the copyright info is already there.
%   addcopyright('My Organization','src')
%
%   % Supply the organization name starting from the "src" folder (found
%   % in the current folder), and use the supplied year but don't
%   % overwrite it if the copyright info is already there.
%   addcopyright('My Organization','private',2001)
%   addcopyright('My Organization','src','2001-2003')
%
%   % Supply the organization name starting from the "src" folder (found
%   % in the current folder), and use the supplied year.  If the
%   % copyright info is already in the M-file, replace it with the new
%   % copyright line.
%   % overwrite it if the copyright info is already there.
%   addcopyright('My Organization','src',2004,true)
%
%   % Replace all instances of the copyright info or add it if it's not
%   % there.  Use the supplied copyright info.
%   addcopyright([],[],[],true, ...
%        'Copyright 1984-2003 MathWorks, makers of MATLAB(R)')

%   Copyright 2003 The MathWorks, Inc.
%   $Revision: 1.2 $ $Date: 2003/09/11 21:25:11 $

%   Author: Raymond Norris (rayn@mathworks.com)

% Parse input arguments
error(nargchk(1,5,nargin))

if nargin<4
    update = false;
end
if nargin<3
    copyright_year = num2str(datestr(date,10));
end
if nargin<2
    folder = pwd;
end

if isempty(folder)
    folder = pwd;
end

if isempty(copyright_year)
    copyright_year = datestr(date,10);
end

if isempty(update)
    update = false;
end

if isnumeric(copyright_year)
    copyright_year = num2str(copyright_year);
end

if nargin<5
    copyright_string = ['%   Copyright ' copyright_year ' ' organization];
else
    copyright_string = ['%   ' copyright_string];
end

dlist = dir(folder);

% Order the listing so that all of the files are first, then the subfolders.
idx = [dlist.isdir];
newdlist = [dlist(~idx); dlist(idx)];

for didx = 1:length(newdlist)
    token_name = newdlist(didx);
    if strcmp(token_name.name,'.') || strcmp(token_name.name,'..') || strcmp(token_name.name,'CVS')
        continue
    end

    if token_name.isdir==true
        addcopyright([folder filesep token_name.name],copyright_year);
        continue
    end

    if length(token_name.name)>2 && strcmp(token_name.name(end-1:end),'.m')
        lAddCopyrightInfo([folder filesep token_name.name],copyright_string,update)
    end
end

function lAddCopyrightInfo(fname,copyright_string,update)

fprintf(1,'Modifying file: %s...',fname)

fid_1 = fopen(fname,'r+');
if fid_1<3
    disp(['Failed to open ' fname ' for modifying.'])
    return
end

[dummy,file,ext] = fileparts(fname);
new_fname = ['new_' file ext];
fid_2 = fopen(new_fname,'w');
if fid_2<3
    disp(['Failed to open ' new_fname ' for writing to.'])
    fclose(fid_1)
    return
end

next_line   = [];
found_token = false;
token_already_embedded = false;
while isempty(next_line) | next_line~=-1
    next_line = fgetl(fid_1);
    if next_line==-1
        break
    elseif strncmp(next_line,'%   Copyright',13)
        token_already_embedded = true;
        found_token = true;
        if update==true
            fprintf(fid_2,'%s\n',copyright_string);
            continue
        end
    end

    if found_token==false
        if strncmp(next_line,'%   $Revision',13)
            fprintf(fid_2,'%s\n',copyright_string);
            found_token = true;
        end
    end
    fprintf(fid_2,'%s\n',next_line);
end

fclose(fid_1);
fclose(fid_2);

if found_token==false
    disp('Could not find Revision token.  Copyright info will NOT be added to the file.')
    delete(new_fname)
elseif token_already_embedded==true && update==false
    disp('Copyright info was found in the file.')
    delete(new_fname)
else
    disp('Done.')
    movefile(new_fname,fname)
end

author: Sable

repo: mcbench-benchmarks


text:
 function [kp,gz,kz,kf,t] = compute_gz_spsp(kp)
%function [kp,gz,kz,kf,t] = compute_gz_spsp(kp)
%Function that computes z gradient waveform accompanying the SPSP pulse,
%using parameters specified in structure kp, generated by kparameterSPSP.m.
%In units of g/cm. Currently can only generate waveforms made of trapezoids.
%Outputs:
%kp: updated trajectory parameter structure
%gz: z gradient waveform in g/cm
%kz,kf: together they specify the trajectory in SPSP k-space (kf-kz space)
%
%Chun-yu Yip, 4/1/2009

%Physical parameters
gam = 26751;                            %rad/sec/g; gyromagnetic ratio
gambar = gam/2/pi;                      %Hz/g

%Trapezoid area calculation
if (kp.dgdtmax*(kp.T/2)/2 <=kp.gmax)             %the lobe is triangular.
    gzarea = kp.T/2 * kp.dgdtmax*(kp.T/2)/2 /2;  %g s /cm
else                                             %this lobe is a trapezoid.
    r = kp.gmax/kp.dgdtmax;                      %s; rise time of gradient
                                                 %to plateau
    gzarea = kp.gmax*r + kp.gmax*(kp.T/2 - 2*r); %g s /cm
end

%Trapezoid generation
gztrap_whole = dotrap(gzarea,kp.gmax,kp.dgdtmax,kp.pointtime);
gztrap_half = dotrap(gzarea/2,kp.gmax,kp.dgdtmax,kp.pointtime);

signz = +1;
gz = [];
for ii = 1:1:kp.Ntraps

  if (ii==kp.Ntraps)
   factor = 1;
  else
   factor = 1;
  end

  signz = -signz;
  gz = [gz factor*signz*gztrap_whole];          %Concatenate trapezoids
                                                %with alternating signs.
end

signz = -signz;
gz = [gz signz*gztrap_half];                    %Add refocusing lobe
gz = gz.';

t = kp.pointtime*[0:1:length(gz)-1].';          %s; time vector

%Do "backward integral" to obtain excitation k-space trajectory
intgz = cumtrapz(gz);
kz = -gambar * kp.pointtime *  (ones(length(gz),1)*intgz(end) - intgz); %/cm
kf = t-t(end);
kp.npnts = length(gz);              %pulse length (number of samples)
kp.pw = kp.pointtime * kp.npnts;    %s; pulse length (actual duration)

if 0 % Display gradient waveform
	figure
	plot(t,gz)
	xlabel('time (sec)')
	ylabel('g/cm')
	title('Gz waveform')
	grid
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function r = DH_bccr(HazeImg,conf)

%  Gaofeng MENG:
%  National Laboratory of Pattern Recognition,
%  Institute of Automation, Academy of Sciences, Beijing 100190
%  Comments and bug reports are welcome.  Email to gfmeng@nlpr.ia.ac.cn
%
%  WORK SETTING:
%  This code has been compiled and tested by using MATLAB R2009a
%
%  For more detials, please see our paper:
%  Gaofeng MENG, Ying WANG, Jiangyong DUAN, Shiming XIANG, Chunhong PAN.
%  Efficient Image Dehazing with Boundary Constraint and Contextual Regularization,
%  ICCV, Sydney, Australia, pp.617-624, 3-6 Dec., 2013.
%
%  Last Modified: Feb. 14, 2014, By Gaofeng MENG

%% param
bccr_method = conf.bccr_method;
bccr_air_wsz = conf.bccr_air_wsz; % airlight window size
bccr_Bou_wsz = conf.bccr_Bou_wsz; % Boundcon window size
bccr_lambda = conf.bccr_lambda; % regularization parameter, the more this parameter, the closer to the original patch-wise transmission

%%
A = bccr_Airlight(HazeImg, bccr_method, bccr_air_wsz);

% calculating boundary constraints
ts = bccr_Boundcon(HazeImg, A, 30, 300, bccr_Bou_wsz);

% refining the estimation of transmission
t = bccr_CalTransmission(HazeImg, ts, bccr_lambda, 0.5); % using contextual information

% dehazing
r = bccr_Dehazefun(HazeImg, t, A, 0.85);

author: AomanHao

repo: Matlab-Image-Dehaze-Enhance


text:
function [Gamma,Xi,L] = fb_Gamma_inference_ehmm(XX,ehmm,residuals,Gamma,Xi,T)
% Inference using forward backward propagation for parallel ehmm chains

% Note that this is an approximation, so there can be small increments of
% the free energy in this step. This is because of two reasons:
% 1. The likelihood for the free energy and for the FB equations is
% computed slightly differently: the term NormWishtrace (related to the
% W covariance) is computed across states for the FB equations, and per
% state in the free energy. That's due to the fact that the FB equations
% are approximate here.
% 2. The weights of how each chain contributes to each time point are *not*
% a linear combination of the Gamma state probabilities, because the way the weight
% for the baseline state is computed is nonlinear (and therefore everything
% is not linear). In fact, the likelihood is also a nonlinear function of
% the state probabilities. For example, we could have the highest
% likelihood for gamma_kt = 0.3, but the forward-backward equations assume
% a linear progression between 0 and 1. A potential solution would be to
% bin the Gammas, and have L states per chain: 0% active, 20% active, etc.,
% but that complicates the definition of the trans prob matrices and their
% priors.
%
% Author: Diego Vidaurre, University of Oxford / Aarhus University (2022)

order = ehmm.train.maxorder; K = ehmm.K; show_warn = true;

if nargin<6, T = size(residuals,1) + order; ttrial = T(1); N = 1;
else, ttrial = T(1); N = length(T);
end % normally this is called for a continuous segment
% if nargin<5 || isempty(Xi), Xi = zeros(ttrial-1-order,K,4); check_improv = false;
% else, check_improv = true;
% end
check_improv = false;

if ehmm.train.acrosstrial_constrained
    Gamma = squeeze(mean(reshape(Gamma,[ttrial-order,N,K]),2));
end

if check_improv
    n_rejected = 0;
    if strcmpi(ehmm.train.stopcriterion,'LogLik')
        c = sum(evalfreeenergy_ehmm(T,Gamma,Xi,ehmm,residuals,XX,[0 1 1 0 0]));
    else
        c = sum(evalfreeenergy_ehmm(T,Gamma,Xi,ehmm,residuals,XX,[1 1 1 0 0]));
    end
    %disp(num2str(loglik))
end

if isempty(Xi), Xi = zeros(ttrial-order-1,K,4); end

for k = randperm(K)

    Lk = zeros(ttrial * N, 2); % N > 1 only if we are averaging across trials

    for j = 1:N % N is 1 except when acrosstrial_constrained
        ind = (1:ttrial) + ttrial * (j-1);
        ind2 = (1:(ttrial-order)) + (ttrial-order) * (j-1);
        Lk(ind,:) = obslike_ehmm(ehmm,Gamma,residuals(ind2,:),XX(ind2,:),k);
    end

    if ehmm.train.acrosstrial_constrained
       Lk = squeeze(sum(reshape(Lk,[ttrial,N,2]),2));
    end

    scale = zeros(ttrial,1);
    alpha = zeros(ttrial,2);
    beta = zeros(ttrial,2);

    alpha(1+order,:) = ehmm.state(k).Pi.*Lk(1+order,:);
    scale(1+order) = sum(alpha(1+order,:));
    alpha(1+order,:) = alpha(1+order,:)/scale(1+order);
    for i = 2+order:ttrial
        alpha(i,:) = (alpha(i-1,:)*ehmm.state(k).P).*Lk(i,:);
        scale(i) = sum(alpha(i,:));		% P(X_i | X_1 ... X_{i-1})
        alpha(i,:) = alpha(i,:)/scale(i);
    end

    scale(scale<realmin) = realmin;

    beta(ttrial,:) = ones(1,2)/scale(ttrial);
    for i = ttrial-1:-1:1+order
        beta(i,:) = (beta(i+1,:).*Lk(i+1,:))*(ehmm.state(k).P')/scale(i);
        beta(i,beta(i,:)>realmax) = realmax;
    end

    if any(isnan(beta(:)))
        warning(['State time courses became NaN (out of precision). ' ...
            'There are probably extreme events in the data. ' ...
            'Using an approximation..' ])
        Gammak = alpha; out_precision = true;
    else
        Gammak = (alpha.*beta); out_precision = false;
    end

    Gammak = rdiv(Gammak,sum(Gammak,2));
    if check_improv, Gammak0 = Gamma(:,k); Xi0k = Xi(:,k,:); end
    Gamma(:,k) = Gammak(order+1:ttrial,1);

    if isempty(Xi), Xi = zeros(ttrial-order-1,K,4); end
    if out_precision
        xi = approximateXi(Gammak(order+1:ttrial,:),size(Gammak,1),ehmm);
        Xi(:,k,:) = reshape(xi,[size(xi,1) size(xi,2)*size(xi,3)]);
        if show_warn
            show_warn = false;
            warning('out of precision when computing xi')
        end
    else
        for i = 1+order:ttrial-1
            t = ehmm.state(k).P .* ( alpha(i,:)' * (beta(i+1,:).*Lk(i+1,:)));
            Xi(i-order,k,:) = t(:)'/sum(t(:));
        end
    end

    if check_improv
        if strcmpi(ehmm.train.stopcriterion,'LogLik')
            c2 = sum(evalfreeenergy_ehmm(T,Gamma,reshape(Xi,T-order-1,K,2,2),...
                ehmm,residuals,XX,[0 1 1 0 0]));
        else
            c2 = sum(evalfreeenergy_ehmm(T,Gamma,reshape(Xi,T-order-1,K,2,2),...
                ehmm,residuals,XX,[1 1 1 0 0]));
        end
        %disp(num2str(loglik2))
        if c2 > c % decrement of fitness, undo
            Gamma(:,k) = Gammak0; Xi(:,k,:) = Xi0k;
            n_rejected = n_rejected + 1;
        else
            c = c2;
        end
    end

end

if check_improv, disp(['%Change rejected: ' num2str(n_rejected/K) ]); end

L = scale; % only the last one

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [thrust, torque] = computePropOpPoint(RPM, rho, d_prop, C_t, C_q)
% [thrust, torque] = computePropOpPoint(RPM, rho, d_prop, C_t, C_q)
%
% Calculates thrust and torque for a series of propeller RPM inputs.
%   Note: Other parameters are considered constant
%   See: https://web.mit.edu/16.unified/www/FALL/thermodynamics/notes/node86.html
%
%
% INPUTS:
%   RPM = [1 x n] (RPM) propeller RPM
%   rho = [scalar] (kg/m^3) air density
%   d_prop = [scalar] (m) propeller diameter
%   C_t = [scalar] () propeller thrust coefficient.
%   C_q = [scalar] () propeller torque coefficient.
%
% OUTPUTS:
%   thrust = [scalar] (N) propeller thrust
%   torque = [scalar] (Nm) propeller torque
%
% Written by Conrad McGreal 2020-1-25

revs = RPM/60 ; % convert to revolution per second

% Calculate thrust
thrust = C_t*rho*revs.^2*d_prop^4 ;

% Calculate torque
torque = C_q*rho*revs.^2*d_prop^5 ;

author: MatthewPeterKelly

repo: OptimTraj


text:
% Init Script
% Remarks:
%   1.  If there is a predefined variable called 'subStreamNumberDefault'
%       it will be used to set the SubStream of the Random Number
%       Generator.
% TODO:
% 	1.  A
% Release Notes
% - 1.0.005     02/08/2019  Royi Avital
%   *   Added Octave Compatibility.
% - 1.0.004     05/07/2018  Royi Avital
%   *   Added case for 'subStreamNumberDefault = 0'.
% - 1.0.003     29/06/2018  Royi Avital
%   *   Clearing all variables but 'subStreamNumberDefault'.
%   *   Setting 'subStreamNumber' only if 'subStreamNumberDefault' not
%       defined.
% - 1.0.002     19/04/2018  Royi Avital
%   *   Added 'markerSizeMedium' (Matches 'markerSizeNormal').
% - 1.0.001     24/11/2016  Royi Avital
%   *   Added figure positions options.
% - 1.0.000     08/11/2016  Royi Avital
%   *   First release.
%
%
%% General Parameters

close('all');
clearvars('-except', 'subStreamNumberDefault');
clc();

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

% See https://stackoverflow.com/questions/2246579
isOctave = (exist('OCTAVE_VERSION', 'builtin') ~= 0);

% set(0, 'DefaultFigureWindowStyle', 'docked');
% defaultLoosInset = get(0, 'DefaultAxesLooseInset');
% set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

figPosSmall     = [100, 100, 0400, 0300];
figPosMedium    = [100, 100, 0800, 0600];
figPosLarge     = [100, 100, 0960, 0720];
figPosXLarge    = [100, 100, 1100, 0825];
figPosX2Large   = [100, 100, 1200, 0900];
figPosX3Large   = [100, 100, 1400, 1050];
figPosDefault   = [100, 100, 0560, 0420];

fontSizeTitle   = 14;
fontSizeAxis    = 12;
fontSizeString  = 12;

lineWidthThin   = 1;
lineWidthNormal = 3;
lineWidthThick  = 4;

markerSizeSmall     = 4;
markerSizeMedium    = 8;
markerSizeNormal    = 8;
markerSizeLarge     = 10;

% https://www.mathworks.com/help/matlab/graphics_transition/why-are-plot-lines-different-colors.html
% https://www.mathworks.com/matlabcentral/answers/160332
mColorOrder = get(groot, 'DefaultAxesColorOrder');

if(isOctave == FALSE)
    randomNumberStream  = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');

    if(exist('subStreamNumberDefault', 'var') && (subStreamNumberDefault ~= 0))
        subStreamNumber = subStreamNumberDefault;
    else
        subStreamNumber = round(sum(clock()));
    end
    set(randomNumberStream, 'Substream', subStreamNumber);
    RandStream.setGlobalStream(randomNumberStream);
end

author: RoyiAvital

repo: StackExchangeCodes


text:
%% SPARSE_FSU uses the fsuClusterMatlab command to run the SPARSE "computation".
%
%  Discussion:
%
%    The arguments to fsuClusterMatlab mean:
%
%    '~/matlab/sparse_parfor' allows us to specify an output directory;
%
%    [] allows us to specify switches to the MOAB queue handler, such as more time.
%
%    'p" indicates this is a "parfor" job;
%
%    'w' indicates that the MATLAB session should wait for completion;
%
%    workers is the number of processors requested;
%
%    @sparse_fun indicates the function we want to execute;
%
%    { n } contains input arguments to the function.
%
%    The function "sparse_fun" must correspond to a MATLAB M-file "sparse_fun.m"
%    and that file must be in MATLAB's path.
%
%    The easiest way to ensure that the file is in MATLAB's path is to
%    create a subdirectory called "matlab" immediately below your home directory,
%    and place the function file there (as well as any other files needed
%    by the job).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2012
%
%  Author:
%
%    John Burkardt
%
  workers = 4;

  tic

  results = fsuClusterMatlab ( '~/matlab/sparse_parfor', [], 'p', 'w', workers, ...
    @sparse_fun, { } );

  toc

  A_sparse = results{1};
%
%  Plot the results.
%
  spy ( A_sparse );

author: johannesgerer

repo: jburkardt-m


text:
function [out] = saturation_2(S,Smax,p1,In)
%saturation_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:  Saturation excess from a store with different degrees of saturation
% Constraints:  1-S/Smax >= 0       prevents numerical issues with complex
%                                   numbers
% @(Inputs):    S    - current storage [mm]
%               Smax - maximum contributing storage [mm]
%               p1   - non-linear scaling parameter [-]
%               In   - incoming flux [mm/d]

% NOTE: When stores are very slightly below or over their maximum, the
% exponent can push this function into regions where no feasible solutions
% exist. The min(max()) combination prevents this from happening.

out = (1- min(1,max(0,(1-S./Smax))).^p1) .*In;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
%SerialLink.JACOB0 Jacobian in world coordinates
%
% J0 = R.jacob0(Q, OPTIONS) is the Jacobian matrix (6xN) for the robot in
% pose Q (1xN), and N is the number of robot joints.  The manipulator
% Jacobian matrix maps joint velocity to end-effector spatial velocity V =
% J0*QD expressed in the world-coordinate frame.
%
% Options::
% 'rpy'     Compute analytical Jacobian with rotation rate in terms of
%           roll-pitch-yaw angles
% 'eul'     Compute analytical Jacobian with rotation rates in terms of
%           Euler angles
% 'trans'   Return translational submatrix of Jacobian
% 'rot'     Return rotational submatrix of Jacobian
%
% Note::
% - The Jacobian is computed in the end-effector frame and transformed to
%   the world frame.
% - The default Jacobian returned is often referred to as the geometric
%   Jacobian, as opposed to the analytical Jacobian.
%
% See also SerialLink.jacobn, jsingu, deltatr, tr2delta, jsingu.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function J0 = jacob0(robot, q, varargin)

    opt.rpy = false;
    opt.eul = false;
    opt.trans = false;
    opt.rot = false;

    opt = tb_optparse(opt, varargin);

	%
	%   dX_tn = Jn dq
	%
	Jn = jacobn(robot, q);	% Jacobian from joint to wrist space

	%
	%  convert to Jacobian in base coordinates
	%
	Tn = fkine(robot, q);	% end-effector transformation
	R = t2r(Tn);
	J0 = [R zeros(3,3); zeros(3,3) R] * Jn;

    if opt.rpy
        rpy = tr2rpy( fkine(robot, q) );
        B = rpy2jac(rpy);
        if rcond(B) < eps
            error('Representational singularity');
        end
        J0 = blkdiag( eye(3,3), inv(B) ) * J0;
    elseif opt.eul
        eul = tr2eul( fkine(robot, q) );
        B = eul2jac(eul);
        if rcond(B) < eps
            error('Representational singularity');
        end
        J0 = blkdiag( eye(3,3), inv(B) ) * J0;
    end

    if opt.trans
        J0 = J0(1:3,:);
    elseif opt.rot
        J0 = J0(4:6,:);
    end

author: Allopart

repo: rbpf-gmapping


text:
% [INPUT]
% y = A vector of floats (-Inf,Inf) of length t representing the residuals.
% a = A float (0,1) representing the significance level (optional, default=0.05).
%
% [OUTPUT]
% h0 = A boolean representing the null hypothesis (residuals are normally distributed) rejection outcome.
% pval = A float [0,1] representing the test p-value.
% stat = A float (-Inf,Inf) representing the normalized test statistic.
%
% [NOTES]
% When residuals are leptokurtic, the Shapiro-Francia test is performed.
% When residuals are platykurtic, the Shapiro-Wilk test is performed.

function [h0,pval,stat] = shapiro_test(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('y',@(x)validateattributes(x,{'double'},{'real' 'finite' 'vector' 'nonempty'}));
        ip.addOptional('a',0.05,@(x)validateattributes(x,{'double'},{'real' 'finite' '>' 0 '<' 1 'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    y = validate_input(ipr.y);
    a = ipr.a;

    nargoutchk(2,3);

    [h0,pval,stat] = shapiro_test_internal(y,a);

end

function [h0,pval,stat] = shapiro_test_internal(x,alpha)

    x = sort(x);

    if (kurtosis(x) > 3)
        [pval,stat] = shapiro_francia(x);
    else
        [pval,stat] = shapiro_wilks(x);
    end

    h0  = (alpha >= pval);

end

function [pval,stat] = shapiro_francia(x)

    n = numel(x);

    x0 = x - mean(x);
    m = norminv(((1:n).' - 0.375) ./ (n + 0.25));
    w = (1 / sqrt(m.' * m)) .* m;
    k = (w.' * x)^2 / (x0.' * x0);

    nu = log(n);
    u1 = log(nu) - nu;
    u2 = log(nu) + (2 / nu);

    mu = -1.27250 + (1.05210 * u1);
    sigma = 1.03080 - (0.26758 * u2);

    stat = (log(1 - k) - mu) / sigma;
    pval = 1 - normcdf(stat,0,1);

end

function [pval,stat] = shapiro_wilks(x)

    n = numel(x);

    x0 = x - mean(x);
    m = norminv(((1:n).' - 0.375) ./ (n + 0.25));
    c = (1 / sqrt(m.' * m)) .* m;
    u = 1 / sqrt(n);

    pc_1 = [-2.7060560  4.434685 -2.071190 -0.147981 0.221157 c(n)];
    pc_2 = [-3.5826330  5.682633 -1.752461 -0.293762 0.042981 c(n - 1)];
    pc_3 = [-0.0006714  0.025054 -0.399780  0.544000];
    pc_4 = [-0.0020322  0.062767 -0.778570  1.382200];
    pc_5 = [ 0.0038915 -0.083751 -0.310820 -1.586100];
    pc_6 = [ 0.0030302 -0.082676 -0.480300];
    pc_7 = [ 0.459 -2.273];

    w = zeros(n,1);
    w(n) = polyval(pc_1,u);
    w(1) = -w(n);

    if (n >= 6)
        off = 3:n-3+1;

        w(n-1) = polyval(pc_2,u);
        w(2) = -w(n-1);

        phi = ((m.' * m) - (2 * m(n)^2) - (2 * m(n-1)^2)) /  (1 - (2 * w(n)^2) - (2 * w(n-1)^2));
    else
        off = 2:n-2+1;

        if (n == 3)
            w(1) = 1 / sqrt(2);
            w(n) = -w(1);

            phi = 1;
        else
            phi = ((m' * m) - (2 * m(n)^2)) / (1 - (2 * w(n)^2));
        end
    end

    w(off) = m(off) ./ sqrt(phi);

    k = (w.' * x)^2 / (x0.' * x0);

    if (n == 3)
        mu = 0;
        sigma = 1;
        kn = 0;
    elseif ((n >= 4) && (n <= 11))
        mu = polyval(pc_3,n);
        sigma = exp(polyval(pc_4,n));
        gamma = polyval(pc_7,n);
        kn = -log(gamma - log(1 - k));
    else
        ln = log(n);
        mu = polyval(pc_5,ln);
        sigma = exp(polyval(pc_6,ln));
        kn = log(1 - k);
    end

    stat = (kn - mu) / sigma;

    if (n == 3)
        pval = (6 / pi()) * (asin(sqrt(k)) - asin(sqrt(0.75)));
    else
        pval = 1 - normcdf(stat,0,1);
    end

end

function y = validate_input(y)

    y = y(:);
    t = numel(y);

    if (t < 5)
        error('The value of ''y'' is invalid. Expected input to be a vector containing at least 5 elements.');
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function c = tone2cent(t1,t2)
% MUSIC.TONE2CENT Returns the number of cents between two semitones.
%    C = MUSIC.TONE2CENT(T) returns the number of cents between T and middle C
%    (C4). T is a semitone interval defined relative to C4.
%
%    C = MUSIC.TONE2CENT(T1,T2) returns the number of cents between T1 and T2.
%    T1 and T2 are musical semitones defined relative to C4. T1 and T2 may be
%    vectors of the same size, or one may be scalar.
%
%    Examples
%       c = music.tone2cent(12);         % returns 1200
%       c = music.tone2cent(-3,-1);      % returns 200
%       c = music.tone2cent([-3 9],-1);  % returns [200 -1000]
%
%    See also music.tone2freq, music.tone2interval, music.tone2note.

%    Author: E. Johnson
%    Copyright 2010 The MathWorks, Inc.

if nargin < 2
    t2 = t1;
    t1 = 0;
end

f1 = music.tone2freq(t1);
f2 = music.tone2freq(t2);

% Could call MUSIC.FREQ2TONE but it is easy to inline this simple calclation.
c  = 1200 * log2( f2 ./ f1 );


text:
function [FRF,FBB] = Receiver(Fopt,NRF)

% randomly generate FRF
[Nt,Ns] = size(Fopt);
FRF = [];
for i = 1:NRF
    FRF = blkdiag(FRF, exp(sqrt(-1) * unifrnd (0,2*pi,[Nt/NRF,1])));
end
FRF = 1/sqrt(Nt)*FRF;

y = [];
while(isempty(y) || abs(y(1)-y(2))>1e-3)
    % fix FRF, optimize FBB
    FBB = pinv(FRF) * Fopt;

    y(1) = norm(Fopt-FRF*FBB,'fro')^2;

    % fix FBB, optimize FRF
    for i = 1:Nt
        m = ceil(i*NRF/Nt);
        FRF(i,m) = 1/sqrt(Nt) * exp( sqrt(-1) * angle( Fopt(i,:)*FBB(m,:)' ) );
    end

    y(2) = norm(Fopt-FRF*FBB,'fro')^2;
end

end

author: yuxianghao


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

function y = StuCVaR4(x)
% General method to compute CVaR for Student distribution with 4 degrees
% of freedom
    f = @(t,n) -(n^(n-2)*(n+t.^2).^(.5-.5*n)*gamma(.5*(n-1))/(2*sqrt(pi)*gamma(.5*n)));
    arg = InverseCDF4(x);
    y = 1 / sqrt(2)*f(arg,4)./x;
end


text:
function value = p24_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P24_F evaluates the integrand for problem 24.
%
%  Dimension:
%
%    DIM_NUM is arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integral Parameters:
%
%    C(1:DIM_NUM) defaults to 0.0.
%
%    A typical, more difficult problem, has
%      C(I) = I**(1/3)
%
%    Call P24_R8VEC to get or set C.
%
%  Integrand:
%
%    F(X) = product (   ( abs ( 4 * X(1:DIM_NUM) - 2 ) + C(1:DIM_NUM) )
%                     / ( 1 + C(1:DIM_NUM) )
%                   )
%
%  Exact Integral:
%
%    1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Joe, Frances Kuo,
%    Remark on Algorithm 659:
%    Implementing Sobol's Quasirandom Seqence Generator,
%    ACM Transactions on Mathematical Software,
%    Volume 29, Number 1, March 2003, pages 49-57.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  c = [];
  c = p24_r8vec ( 'G', 'C', dim_num, c );

  value(1:point_num) = 0.0;

  for point = 1 : point_num

    value(point) = prod ...
    ( ...
        ( abs ( 4.0 * x(1:dim_num,point) - 2.0 ) + c(1:dim_num)' ) ...
      ./ ( 1.0 + c(1:dim_num)' ) ...
    );
  end

  p24_i4 ( 'I', '#', point_num );

  return
end


text:
function x= cauchyinv(p, varargin)

% USAGE:       x= cauchyinv(p, a, b)
%
% Inverse of the Cauchy cumulative distribution function (cdf), x= a + b*tan(pi*(p-0.5)).
%
% ARGUMENTS:
% p (0<=p<=1) might be of any dimension.
% a (default value: 0.0) must be scalars or size(p).
% b (b>0, default value: 1.0) must be scalars or size(p).
%
% EXAMPLE:
% p= 0:0.01:1;
% plot(cauchyinv(p), p);
%
% SEE ALSO:    cauchycdf, cauchyfit, cauchypdf, cauchyrnd.
%
% Copyright (C) Peder Axensten <peder at axensten dot se>
%
% HISTORY:
% Version 1.0, 2006-07-10.
% Version 1.1, 2006-07-26.
% - Added cauchyfit to the cauchy package.
% Version 1.2, 2006-07-31:
% - cauchyinv(0, ...) returned a large negative number but should be -Inf.
% - Size comparison in argument check didn't work.
% - Various other improvements to check list.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	% Default values
	a=	0.0;
	b=	1.0;

	% Check the arguments
	if(nargin >= 2)
		a=	varargin{1};
		if(nargin == 3)
			b=			varargin{2};
			b(b <= 0)=	NaN;	% Make NaN of out of range values.
		end
	end
	if((nargin < 1) || (nargin > 3))
		error('At least one argument, at most three!');
	end

	p(p < 0 | 1 < p)=	NaN;

	% Calculate
	x=			a + b.*tan(pi*(p-0.5));

	% Extreme values.
	if(numel(p) == 1), 	p= repmat(p, size(x));		end
	x(p == 0)=	-Inf;
	x(p == 1)=	Inf;
end


text:
function Gamma = vorticity(L)
% vorticity Gamma

Gamma = L.rotationRate ./ L.strainRate;

author: mtex-toolbox

repo: mtex

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


text:
function [ p, seed ] = triangle_sample ( t, n, seed )

%*****************************************************************************80
%
%% TRIANGLE_SAMPLE returns random points in a triangle.
%
%  Discussion:
%
%    In order to compare the results of this program with the
%    C++ and FORTRAN90 versions, it is useful to use R8VEC_UNIFORM_01
%    to generate the random numbers.  However, it is much faster
%    to rely on MATLAB's RAND routine to do this.  So if speed
%    is a consideration, replace the calls to R8VEC_UNIFORM_01
%    as indicated.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Input, integer N, the number of points to generate.
%
%    Input, integer SEED, a seed for the random number generator.
%    Only needed if R8VEC_UNIFORM_01 is used.
%
%    Output, real P(2,N), random points in the triangle.
%
%    Output, integer SEED, a seed for the random number generator.
%
  dim_num = 2;
%
%  For fast execution, call RAND:
%
% alpha = rand ( 1, n );
%
%  For comparison with  F90 and C++, call R8VEC_UNIFORM:
%
  [ alpha, seed ] = r8vec_uniform_01 ( n, seed );
%
%  Interpret R as a percentage of the triangle's area.
%
%  Imagine a line L, parallel to side 1, so that the area between
%  vertex 1 and line L is R percent of the full triangle's area.
%
%  The line L will intersect sides 2 and 3 at a fraction
%  ALPHA = SQRT ( R ) of the distance from vertex 1 to vertices 2 and 3.
%
  alpha(1:n) = sqrt ( alpha(1:n) );
%
%  Determine the coordinates of the points on sides 2 and 3 intersected
%  by line L.
%
  for dim = 1 : dim_num
    p12(dim,1:n) = ( 1.0 - alpha(1:n) ) * t(dim,1) ...
                         + alpha(1:n)   * t(dim,2);

    p13(dim,1:n) = ( 1.0 - alpha(1:n) ) * t(dim,1) ...
                         + alpha(1:n)   * t(dim,3);
  end
%
%  Now choose, uniformly at random, a point on the line L.
%
%  For fast execution, call RAND:
%
% alpha = rand ( 1, n );
%
%  For comparison with  F90 and C++, call R8VEC_UNIFORM:
%
  [ alpha, seed ] = r8vec_uniform_01 ( n, seed );

  for dim = 1 : dim_num
    p(dim,1:n) = ( 1.0 - alpha(1:n) ) .* p12(dim,1:n) ...
                       + alpha(1:n)   .* p13(dim,1:n);
  end

 return
end


text:
function [str] = printDVManeuversMFMSToTextbox(hDvManInfoText,waypoints,dVDepartVectNTW,deltaVVectNTW,eOrbit,orbitsIn,paddLen,form)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

    hRule = getHRule();
    str = {};

    sumDv = 0;
    for(i=1:length(waypoints)-1) %#ok<*NO4LP>
        if(i==1)
            dvVect = dVDepartVectNTW;
            ta = eOrbit(8);
        else
            dvVect = deltaVVectNTW(:,i-1);
            ta = orbitsIn(i-1,6);
        end

        str{end+1} = ['Burn Information to Depart ', cap1stLetter(waypoints{i}.name)];
        str{end+1} = hRule;
        str{end+1} = [paddStr('Total Delta-V = ',paddLen), num2str(norm(dvVect), form), ' km/s'];
        str{end+1} = [paddStr('Prograde Delta-V = ',paddLen), num2str(1000*dvVect(1), form), ' m/s'];
        str{end+1} = [paddStr('Orbit Normal Delta-V = ',paddLen), num2str(1000*dvVect(2), form), ' m/s'];
        str{end+1} = [paddStr('Radial Delta-V = ',paddLen), num2str(1000*dvVect(3), form), ' m/s'];
        str{end+1} = '---------------------';
        str{end+1} = [paddStr('Burn True Anomaly = ',paddLen), num2str(rad2deg(ta), form), ' deg'];
        str{end+1} = hRule;

        sumDv = sumDv + norm(dvVect);
    end

    str{end+1} = [paddStr('Total Mission Delta-V = ',paddLen), num2str(sumDv, form), ' km/s'];

    set(hDvManInfoText,'String',str);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [r,spherCent]=osculatingSpher4LatLon(latLon,a,f,is2D)
%%OSCULATINGSPHER4LATLON Given the latitude and longitude of a point on an
%           ellipsoid, determine the radius and center location of an
%           osculating sphere going through that point. The osculating
%           sphere is such that the azimuth equals the longitude and the
%           elevation equals the latitude of the point on the ellipsoid.
%           Additionally, the local tangent planes at the point are the
%           same. The radius of curvature of the sphere at that point
%           should nominally equal the radius of curvature of the ellipsoid
%           at that point. However, since the ellipsoid has two radii of
%           curvature, the geometric mean is used (the mean radius of
%           curvature). On the other hand, if is2D is true, then just the
%           radius of curvature in the meridian will be used, which is what
%           perfectly matches a cut of the ellipsoid through the origin and
%           pole (A 2D ellipse in the x-z plane).
%
%INPUTS: latLon The 2X1 point of the format [latitude;longitude] in radians
%               where the osculating sphere should touch the ellipsoid.
%             a The semi-major axis length of the ellipsoid. If this
%               argument is omitted or an empty matrix is passed, the value
%               in Constants.WGS84SemiMajorAxis is used.
%             f The flattening factor of the ellipsoid. If this argument is
%               omitted or an empty matrix is passed, the value in
%               Constants.WGS84Flattening is used.
%          is2D Indicates whether one just cares about a 2D cut of the
%               ellipsoid. If so, then the radius used will be the radius
%               of curvature in the meridian. The default if omitted or an
%               empty matrix is passed is false.
%
%OUTPUTS: r The radius of the osculating sphere.
% spherCent The location of the center of the osculating sphere with
%           respect to the center of the ellipsoid.
%
%The formulae in Section 3 of [1] are used.
%
%EXAMPLE:
%Here, we consider a 2D slice of an ellipsoid that goes through
%longitude=0, we choose a point. We then find the osculating sphere. The
%ellipse cut of the ellipsoid and the circle cut of the sphere are then
%plotted with the point. We ALSO plot the circle cut obtained with
%is2D=true. Since we are only considering a 2D cut, one will see that that
%circle is actually better aligned with the curvature of the ellipse at
%that point. However, if one were to take a 2D cut in other directions, it
%would appear to be a worse fit.
% a=20;%Semi-major axis.
% f=0.5;%Flattening factor.
% b=a*(1-f);%The semi-minor axis of the ellipsoid.
%
% figure(1)
% clf
% hold on
% %Draw a 2D cut of the ellipsoid.
% A=inv([a^2,0;
%        0,b^2]);
% drawEllipse([0;0],A,1,'linewidth',4)
% latLon=[50;0]*(pi/180);
% oscPt=ellips2Cart([latLon;0],a,f);
% [r,spherCent]=osculatingSpher4LatLon(latLon,a,f,false);
% [r2D,spherCent2D]=osculatingSpher4LatLon(latLon,a,f,true);%is2D=true
% %Draw the both 2D cuts of the sphere:
% A=inv([r^2,0;
%        0,r^2]);
% drawEllipse([spherCent(1);spherCent(3)],A,1,'-.','linewidth',4)
% A=inv([r2D^2,0;
%        0,r2D^2]);
% drawEllipse([spherCent2D(1);spherCent2D(3)],A,1,'--','linewidth',2)
%
% %Note that spherCent(2)==0.
% scatter(oscPt(1),oscPt(3),100,'filled')
% legend('Ellipsoid Cut','Osculating Sphere','Osculating Sphere 2D','location','southwest')
% h1=xlabel('x');
% h2=ylabel('y');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] P. Williams and D. Last, "On Loran-C time-difference to co-ordinate
%    converters," in Proceedings of the 32nd Annual Convention & Technical
%    Symposium of the International Loran Association, Boulder, CO, 3-7
%    Nov. 2003.
%
%August 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

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

if(nargin<3||isempty(f))
    f=Constants.WGS84Flattening;
end

if(nargin<2||isempty(a))
    a=Constants.WGS84SemiMajorAxis;
end

%The squared eccentricity of the ellipsoid.
e2=f*(2-f);

lat=latLon(1);%Latitude.
lon=latLon(2);

Phi=lat;

denomTerm=sqrt(1-e2*sin(Phi)^2);
%Equation 3.17 in [1]. The radius of curvature in the meridian.
M0=a*(1-e2)/denomTerm^3;

if(is2D)
    r=M0;
else
    %Equation 3.18 in [1]. The radius of curvature in the prime vertical.
    N0=a/denomTerm;

    %Equation 3.16 in [1].
    r=sqrt(M0*N0);
end

xCartEllips=ellips2Cart([latLon(1:2);0],a,f);
xCartSpher=spher2Cart([r;lon;lat]);

spherCent=-(xCartSpher-xCartEllips);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ x, w ] = ncc_nested_rule ( n, x_min, x_max )

%*****************************************************************************80
%
%% NCC_NESTED_RULE computes the nested Newton Cotes Closed rule.
%
%  Discussion:
%
%    Nested version of Newton Cotes Closed rule on [-1,+1].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), W(N), the points and weights.
%

%
%  Get N equally spaced points that, in general, include [ 0, 1];
%
  if ( n == 1 )
    x = 0.5;
  elseif ( n == 2 )
    x = [ 0.5, 0.0 ];
  elseif ( n == 3 )
    x = [ 0.5, 0.0, 1.0 ];
  else
    x = vdc_sequence ( n - 2 );
    x(4:n) = x(2:n-2);
    x(2) = 0.0;
    x(3) = 1.0;
  end
%
%  Linear transformation from [0,1] to [-1,+1].
%
  x = 2.0 * x - 1.0;
%
%  Compute weights.
%
  x_min = -1.0;
  x_max = +1.0;

  w = nc_compute ( n, x_min, x_max, x );

  return
end


text:
function drawEllipse(h, x, P, c, ns, NP)

% DRAWELLIPSE Draw 2D ellipse or 3D ellipsoid.
%   DRAWELLIPSE(H,X,P) redraws ellipse or ellipsoid in handle H with mean X
%   and covariance P.
%
%   DRAWELLIPSE(...,C) allows the specification of the ellipse color C. If
%   C is not given or is empty, the color is unchanged.
%
%   DRAWELLIPSE(...,C,NS,NP) allows entering the number of sigmas NS and
%   the number of points NP. The default values are NS=3 and NP=10. Use
%   C=[] if you do not want to modify the ellipse color.

%   Copyright 2009 Joan Sola @ LAAS-CNRS.

if nargin < 6,         NP = 10;
    if nargin < 5,     ns = 3;
    end
end

if numel(x) == 2

    [X,Y] = cov2elli(x,P,ns,NP);
    set(h,'xdata',X,'ydata',Y,'vis','on')

elseif numel(x) == 3

    [X,Y,Z] = cov3elli(x,P,ns,NP);
    set(h,'xdata',X,'ydata',Y,'zdata',Z,'vis','on')

else
    error('??? Size of vector ''x'' not correct.')
end

if nargin > 3 && ~isempty(c)
    set(h,'color',c)
end

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function writeLogData(fid,data,format)

% WRITELOGDATA  Write log data in file.
%   WRITELOGDATA(FID,DATA) writes the line of data DATA in hte log file
%   identified by FID. FID is created by the function OPENLOGFILE.
%
%   See also OPENLOGFILE, WRITELOGHEADER, CLOSELOGFILE, FPRINTF.

if nargin == 2
    n = length(data);
    format = '%d';
    for i = 1:n-1
        format = [format ' %f'];
    end
    format = [format '\n'];
end

fprintf(fid,format,data);

end
% return

%%
% fid = fopen('lala.log','w');
% writeLogHeader(fid,'f x');
% writeLogData(fid,3);
% fclose(fid);

%%

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE


text:
function mtrTrackVis2PDB(volFile, trkFile, pdbFile)

% File format conversion from the .trk file format (TrackVis) to pdb files
%
% Parameters
% ----------
% volFile: A full-path to a nifti file, from which the qform will be read
% trkFile: The trk file to be converted.
% pdbFile: The output.
%
% Returns
% -------
% Nothing.

% Load vol file to get AcPc xform
vol = niftiRead(volFile);
xformTo = vol.qto_xyz;
xformFrom = abs(vol.qto_ijk);
xformFrom(1:3,4)=0;
xformToAcPc = xformTo*xformFrom;

% Get mm dimensions of image
dim_mm = vol.dim(:).*abs(diag(vol.qto_xyz(1:3,1:3)));
clear vol;

% Load trks file
fg = read_trk_to_fg(trkFile);

% Remove really short fibers (<10 points)
%fg = dtiCleanFibers(fg, [], [], 11)

fg = dtiXformFiberCoords(fg,xformToAcPc);
fg = dtiClearQuenchStats(fg);
fg = dtiCreateQuenchStats(fg,'Length','Length', 1);
mtrExportFibers(fg,pdbFile,eye(4));

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = createMengerSponge()
%CREATEMENGERSPONGE Create a cube with an inside cross removed
%
%   [n, e, f] = createMengerSponge;
%   Main use is to test possibility of drawing polyhedra with complex faces
%   (polygonal faces with holes)
%
%   Example
%   [n, e, f] = createMengerSponge;
%   drawMesh(n, f);
%
%   See also
%   meshes3d, drawMesh
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2007-10-18
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%   HISTORY
%   2008-10-17 finishes implementation

nodes =[...
    ... % main cube corners (1->8)
    0 0 0; ...
    3 0 0; ...
    0 3 0; ...
    3 3 0; ...
    0 0 3; ...
    3 0 3; ...
    0 3 3; ...
    3 3 3; ...
    ... % outer cube inner face corners
    1 1 0; ... % face z=0 (9->12)
    2 1 0; ...
    1 2 0; ...
    2 2 0; ...
    1 1 3; ... % face z=3 (13->16)
    2 1 3; ...
    1 2 3; ...
    2 2 3; ...
    1 0 1; ... % face y=0 (17->20)
    2 0 1; ...
    1 0 2; ...
    2 0 2; ...
    1 3 1; ... % face y=3 (21->24)
    2 3 1; ...
    1 3 2; ...
    2 3 2; ...
    0 1 1; ... % face x=0 (25->28)
    0 2 1; ...
    0 1 2; ...
    0 2 2; ...
    3 1 1; ... % face x=3 (29->32)
    3 2 1; ...
    3 1 2; ...
    3 2 2; ...
    ... % inner cube corners  (33->40)
    1 1 1; ...
    2 1 1; ...
    1 2 1; ...
    2 2 1; ...
    1 1 2; ...
    2 1 2; ...
    1 2 2; ...
    2 2 2; ...
    ];

edges = [...
    1 2;1 3;2 4;3 4;5 6;5 7;6 8;7 8;1 5;2 6;3 7;4 8;... % outer cube
    9 10;9 11;10 12;11 12;13 14;13 15;14 16;15 16; ...
    17 18;17 19;18 20;19 20; 21 22;21 23;22 24;23 24; ...
    25 26;25 27;26 28;27 28; 29 30;29 31;30 32;31 32; ...
    33 34;33 35;34 36;35 36; 37 38;37 39;38 40;39 40; ... % inner cube
    33 37;34 38;35 39;36 40; ...
     9 33;10 34;11 35;12 36; ... % parallel to xy
    13 37;14 38;15 39;16 40; ...
    17 33;18 34;19 37;20 38; ... % parallel to yz
    21 35;22 36;23 39;24 40; ...
    25 33;26 35;27 37;28 39; ... % parallel to xz
    29 34;30 36;31 38;32 40; ...
    ];

% Alternative definition for faces:
%     [1 2 4 3 NaN 9  11 12 10], ...
%     [5 6 8 7 NaN 13 15 16 14], ...
%     [1 5 7 3 NaN 25 26 28 27], ....
%     [2 6 8 4 NaN 29 30 32 31], ...
%     [1 2 6 5 NaN 17 18 20 19], ...
%     [3 4 8 7 NaN 21 22 24 23], ...

faces = {...
    ... % 6 square faces with a square hole
    [1 2 4 3 1 9  11 12 10  9], ...
    [5 6 8 7 5 13 15 16 14 13], ...
    [1 5 7 3 1 25 26 28 27 25], ....
    [2 6 8 4 2 29 30 32 31 29], ...
    [1 2 6 5 1 17 18 20 19 17], ...
    [3 4 8 7 3 21 22 24 23 21], ...
    ... % faces orthogonal to XY plane, parallel to Oz axis
    [ 9 10 34 33], [ 9 11 35 33], [10 12 36 34], [11 12 36 35], ...
    [13 14 38 37], [13 15 39 37], [14 16 40 38], [15 16 40 39], ...
    ... % faces orthogonal to YZ plane, parallel to Oy axis
    [17 18 34 33], [17 19 37 33], [18 20 38 34], [19 20 38 37], ...
    [21 22 36 35], [21 23 39 35], [22 24 40 36], [23 24 40 39], ...
    ...% faces orthogonal to the YZ plane, parallel to Ox axis
    [25 33 35 26], [25 33 37 27], [26 35 39 28], [27 37 39 28], ...
    [29 30 36 34], [29 31 38 34], [30 32 40 36], [31 32 40 38], ...
    };

% format output
varargout = formatMeshOutput(nargout, nodes, edges, faces);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [yi, ypi, yppi] = cubiconv(x, y, xi)

% CUBICONV Cubic convolution interpolation
%    CUBICONV(X,Y,XI) interpolates to find YI, the value of
%    the underlying function Y at the points in the uniform array
%    XI, using cubic convolution interpolation.  X and Y must be
%    vectors of length N.
%
%    [YI,YPI,YPPI] = CUBICONV() also returns the interpolated
%    quartic derivative and cubic second derivative of the underlying
%    function Y at points XI.

% Joe Henning - Fall 2011

% Cubic Convolution Interpolation for Digital Image Processing
% Robert G. Keys
% IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP-29, No. 6
% December, 1981

n = length(x);
%h = (x(n) - x(1))/n;

for i = 1:length(xi)
   if (xi(i) < x(1) || xi(i) > x(n))
      fprintf('??? Bad x input to cubiconv ==> x(1) <= x <= x(n)\n');
      yi(i) = NaN;
      continue;
   end

   % Find the right place in the table by means of a bisection.
   klo = 1;
   khi = n;
   while (khi-klo > 1)
      k = fix((khi+klo)/2.0);
      if (x(k) > xi(i))
         khi = k;
      else
         klo = k;
      end
   end

   h = x(khi) - x(klo);
   if (h == 0.0)
      fprintf('??? Bad x input to cubiconv ==> x values must be distinct\n');
      yi(i) = NaN;
      ypi(i) = NaN;
      yppi(i) = NaN;
      continue;
   end

   km1 = klo - 1;
   kp1 = khi + 1;
   if (km1 < 1)
      a = 3*y(klo) - 3*y(khi) + y(khi+1);
   else
      a = y(km1);
   end
   if (kp1 > n)
      d = 3*y(khi) - 3*y(klo) + y(klo-1);
   else
      d = y(kp1);
   end
   b = y(klo);
   c = y(khi);

   % Evaluate cubic polynomial
   t = (xi(i) - x(klo))/h;
   t2 = t*t;
   t3 = t2*t;
   h2 = h*h;
   c00 = (-t3 + 2*t2 - t)/2.0;
   c10 = (3*t3 - 5*t2 + 2)/2.0;
   c20 = (-3*t3 + 4*t2 + t)/2.0;
   c30 = (t3 - t2)/2.0;

   yi(i) = a*c00 + b*c10 + c*c20 + d*c30;

   % Differentiate to find the second-order interpolant
   c00 = (-3*t2 + 4*t - 1)/2.0;
   c10 = (9*t2 - 10*t)/2.0;
   c20 = (-9*t2 + 8*t + 1)/2.0;
   c30 = (3*t2 - 2*t)/2.0;

   ypi(i) = a*c00/h + b*c10/h + c*c20/h + d*c30/h;

   % Differentiate to find the first-order interpolant
   c00 = (-6*t + 4)/2.0;
   c10 = (18*t - 10)/2.0;
   c20 = (-18*t + 8)/2.0;
   c30 = (6*t - 2)/2.0;

   yppi(i) = a*c00/h2 + b*c10/h2 + c*c20/h2 + d*c30/h2;
end

author: Sable


text:
% SYNTAX:
% intensity = hmrR_PreprocessIntensity_MedianFilter( intensity )
%
% UI NAME:
% hmrR_PreprocessIntensity_MedianFilter
%
% DESCRIPTION:
% Applies a median filter to data to remove huge spikes.
%
% INPUT:
% intensity - SNIRF data type where the d matrix is intensity
%
% OUTPUT:
% intensity - SNIRF data type where the d matrix is intensity
%
% USAGE OPTIONS:
% Intensity_to_Intensity: d = hmrR_PreprocessIntensity_MedianFilter(data)

function d = hmrR_PreprocessIntensity_MedianFilter( intensity)

for ii=1:length(intensity)
    d = intensity(ii).GetDataTimeSeries();

    for j = 1:size(d,2)
        foo = d(:,j);

        xdata = (1:length( foo))';
        d(:,j) = interp1(xdata(~isnan( foo)), foo(~isnan( foo)),xdata,'spline');

        new_signal = zeros(size(foo));
        new_signal(1) = foo(1);
        new_signal(end) = foo(end);
        for tp = 2:length(foo)-1
            values = foo([tp-1,tp,tp+1]);
            median_value = median(values);
            new_signal(tp) = median_value;
        end
        d(:,j) = new_signal;

    end

    intensity(ii).SetDataTimeSeries(d);
end

author: BUNPC


text:
function [WMM,VADS] = CalMetric(PopObj,W)
% Calculate the weighted min-max metric and the vector angle distance
% scaling metric between each solution and vector

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

	N  = size(PopObj,1);
    NW = size(W,1);

    %% Normalization
    Zmax   = max(PopObj,[],1);
    Zmin   = min(PopObj,[],1);
    PopObj = (PopObj-repmat(Zmin,N,1))./(repmat(Zmax-Zmin,N,1));
    W      = 1./(W-repmat(Zmin,NW,1)+eps);

    %% Calculate WMM
    WMM = zeros(N,NW);
    for i = 1 : NW
        WMM(:,i)  = max(PopObj.*repmat(W(i,:)./norm(W(i,:)),N,1),[],2);
    end

    %% Calculate VADS
    if nargout > 1
        VADS  = zeros(N,NW);
        NormP = sqrt(sum(PopObj.^2,2));
        for i = 1 : NW
            VADS(:,i) = NormP./(PopObj*W(i,:)'./NormP).^100;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function obj = setType(obj, type)
    % Sets the type (either LINEAR or NONLINEAR) of the NLP function
    %
    % Parameters:
    % type: the type of the NLP function @type char

    % specify the function type
    if strcmpi(type, 'Linear')
        obj.Type = 'Linear';
    elseif strcmpi(type, 'Nonlinear')
        obj.Type = 'Nonlinear';
    else
        error('NlpFunction:incorrectFunctionType',...
            'Unspecified function type detected.\n');
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
% peq.m - Parametric EQ with matching gain at Nyquist frequency
% Author: Sophocles J. Orfanidis,
% Usage:  [b, a] = peq(G0, G, w0, Dw, NdB)
%
% G0 = reference gain at DC in dB
% G  = boost/cut gain in dB
% GB = bandwidth gain in dB
%
% w0 = center frequency in Hz
% Dw = bandwidth in Hz
% NdB = amount of db less than the peak/notch gain as the bandwidth frequencies
% b  = [b0, b1, b2] = numerator coefficients
% a  = [1,  a1, a2] = denominator coefficients

% Modified by Arvind using the equations from IIR Filter Design Chapter in
% the book
function [b, a, beta] = peq(G0, G, w0, Dw, NdB)

%Convert from Decibels to real values
GB = G * NdB;
G  = (10^(G/20));
GB = (10^(GB/20));
G0 = (10^(G0/20));

%Convert absolute frequencies to rads/sec
%w0 = w0*2*pi/fs;
%Dw = Dw*2*pi/fs;

beta = sqrt((GB^2-G0^2)/(G^2-GB^2)) * tan(Dw/2);
b = [(G0+G*beta)/(1+beta) -(2*G0*cos(w0))/(1+beta)  (G0-G*beta)/(1+beta)];
a = [1 -(2*cos(w0))/(1+beta) (1-beta)/(1+beta)];

%***********************
% double beta;
%
% beta = sqrt((GB[0]^2-G0[0]^2)/(G^2[0]-GB[0]^2)) * tan(Dw[0]/2)
% Num[0] = (G0[0]+G[0]*beta)/(1+beta);
% Num[1] =  -(2*G0[0]*cos(w0[0]))/(1+beta) ;
% Num[2] = (G0[0]-G[0]*beta)/(1+beta);
%
% Den[0] = 1;
% Den[1] = -(2*cos(w0[0]))/(1+beta);
% Den[2] = (1-beta)/(1+beta);
%*************************

repo: mcbench-benchmarks


text:
%% Setting of the problem
global s
pde = fracLapdata9;
% pde = checkboarddata;
pde.L = 1;
option.maxIt = 4;
option.maxN = 1e6;
option.elemType = 'P1P2';
option.solver = 'mg';
option.gNquadorder = 4;
option.tol = 1e-6;
[node,elem] = squaremesh([-1 1 -1 1],0.5);
bdFlag = setboundary(node,elem,'Dirichlet');

%% s = 0.2
s = 0.2;
femfracLap(node,elem,pde,bdFlag,option);

%% s = 0.4
s = 0.4;
femfracLap(node,elem,pde,bdFlag,option);

%% s = 0.6
s = 0.6;
femfracLap(node,elem,pde,bdFlag,option);

%% s = 0.8
s = 0.8;
femfracLap(node,elem,pde,bdFlag,option);

repo: ifem


text:
function r = uplus(p)
%UMINUS: this function applies a unary minus to the x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

r=p;

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [ a, b ] = p04_lim ( dim_num )

%*****************************************************************************80
%
%% P04_LIM returns the integration limits for problem 04.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end


text:
function M = prepare_mesh_fittemplate(headshape,template)

% PREPARE_MESH_FITTEMPLATE computes an affine transformation matrix between 2 point clouds
%
% This function relies on cpd toolbox from  Myronenko, see https://sites.google.com/site/myronenko/research/cpd
%
%
% See also FT_PREPARE_MESH

% Copyright (C) 2019, Simon Homoelle
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% add toolbox cpd
ft_hastoolbox('cpd', 2);

%
opt.corresp = 0;
opt.method  = 'affine';
opt.max_it = 100;
opt.fgt = 0;
opt.tol = 10e-12;
opt.outliers = 0.0;
opt.outliers = 0;
[transform,~] = cpd_register(headshape,template, opt);

M = eye(4,4);
M(1:3,1:3) = transform.R;
M(1:3,4)   = transform.t;

end

author: spm


text:
function opt_checkTypes(opt, optDefaults, structname)
%OPT_CHECKTYPE - Check types of the values in property/value struct
%
%Synopsis:
%  opt_checkTypes(OPT, PROPSPEC)
%
%Arguments:
%  PROPSPEC: PROPSPECLIST - Property specification list with type definition,
%      i.e., CELL of size [N 3], with the first column all being strings
%      (property names), the second column holding the default values (not
%      used in this function), and the third column specifying the expected
%      type of the values. If a third column does not exist, this function
%      just exists.
%
%Returns:
%  nothing (just throws errors)
%
%Description:
%  This function ensures that the values in the property/value struct OPT
%  comply with the types that are specified in the property specification
%  PROCSPEC. For a list of possible specifications, see misc_checkType.
%
%See also misc_checkType, opt_setDefaults.
%
%Examples:
%  props= {'LineColor', 'k', 'CHAR[1]|DOUBLE[3]';
%          'LineWidth', 0.7, 'DOUBLE';
%          'Colormap', copper(51), 'DOUBLE[- 3]';
%          'Transform', eye(2,2), 'DOUBLE[2 2]'};
%
%  opt= struct('Colormap', hot(5));
%  opt_checkTypes(opt, props)
%
%  opt= struct('LineColor', [1 0 0.5 0.5]);
%  % Should throw an error
%  opt_checkTypes(opt, props)
%
%  opt= struct('LineWidth', 'Very Thick');
%  % Should throw an error
%  opt_checkTypes(opt, props)
%
%  opt= struct('Transform', randn(2,2));
%  opt_checkTypes(opt, props)

% 06-2012 Benjamin Blankertz

if size(optDefaults, 2) < 3,
  % No types specified
  return
end

for k= 1:size(optDefaults, 1),
  fld= optDefaults{k,1};
  if isfield(opt, fld),
    fielddisplayname= [structname '.' fld];
    for n= 1:length(opt),
      [ok, msg]= misc_checkType(opt(n).(fld), optDefaults{k,3}, fielddisplayname);
      if ~ok,
        error(msg);
      end
    end
  end
end

author: bbci