text:
function [n_best P_best V C]=fast_oopsi(F,V,P)
% this function solves the following optimization problem:
% (*) n_best = argmax_{n >= 0} P(n | F)
% which is a MAP estimate for the most likely spike train given the
% fluorescence signal.  given the model:
%
% <latex>
% \begin{align}
% C_t &= \gamma*C_{t-1} + n_t, \qquad & n_t & \sim \text{Poisson}(n_t; \lamda_t \Delta)
% F_t &= \alpha(C_t + \beta) + \sigma \varepsilon_t, &\varepsilon_t &\sim \mathcal{N}(0,1)
% \end{align}
% </latex>
%
% if F_t is a vector, then 'a' is a vector as well
% we approx the Poisson with an Exponential (which means we don't require integer numbers of spikes).
% we take an "interior-point" approach to impose the nonnegative contraint on (*).
% each step is solved in O(T)
% time by utilizing gaussian elimination on the tridiagonal hessian, as
% opposed to the O(T^3) time typically required for non-negative
% deconvolution.
%
% Input---- only F is REQUIRED.  the others are optional
% F:        fluorescence time series (can be a vector (1 x T) or a matrix (Np x T)
%
% V.        structure of algorithm Variables
%   Ncells: # of cells within ROI
%   T:      # of time steps
%   Npixels:# of pixels in ROI
%   dt:     time step size, ie, frame duration, ie, 1/(imaging rate)
%   n:      if true spike train is known, and we are plotting, plot it (only required is est_a==1)
%   h:      height of ROI (assumes square ROI) (# of pixels) (only required if est_a==1 and we are plotting)
%   w:      width of ROI (assumes square ROI) (# of pixels) (only required if est_a==1 and we are plotting)
%
%   THE FOLLOWING FIELDS CORRESPOND TO CHOICES THAT THE USER MAKE
%
%   fast_poiss:     1 if F_t ~ Poisson, 0 if F_t ~ Gaussian
%   fast_nonlin:    1 if F_t is a nonlinear f(C_t), and 0 if F_t is a linear f(C_t)
%   fast_plot:      1 to plot results after each pseudo-EM iteration, 0 otherwise
%   fast_thr:       1 if thresholding inferred spike train before estiamting {a,b}
%   fast_iter_max:  max # of iterations of pseudo-EM  (1 to use default initial parameters)
%   fast_ignore_post: 1 to keep iterating pseudo-EM even if posterior is not increasing, 0 otherwise
%
%   THE BELOW FIELDS INDICATE WHETHER ONE WANTS TO ESTIMATE EACH OF THE
%   PARAMETERS. IF ANY IS SET TO ZERO, THEN WE DO NOT TRY TO UPDATE THE
%   ORIGINAL ESTIMATE, GIVEN EITHER BY THE USER, OR THE INITIAL ESTIMATE
%   FROM THE CODE
%
%   est_sig:    1 to estimate sig
%   est_lam:    1 to estimate lam
%   est_gam:    1 to estimate gam
%   est_b:      1 to estimate b
%   est_a:      1 to estimate a
%
% P.        structure of neuron model Parameters
%
%   a:      spatial filter
%   b:      background fluorescence
%   sig:    standard deviation of observation noise
%   gam:    decayish, ie, tau=dt/(1-gam)
%   lam:    firing rate-ish, ie, expected # of spikes per frame
%
% Output---
% n_best:   inferred spike train
% P_best:   inferred parameter structure
% V:        structure of Variables for algorithm to run

%% initialize algorithm Variables
starttime   = cputime;
siz         = size(F);      if siz(2)==1, F=F'; siz=size(F); end
j=0;
% variables determined by the data
if nargin < 2,              V   = struct;       end
if ~isfield(V,'Ncells'),    V.Ncells = 1;       end     % # of cells in image
if ~isfield(V,'T'),         V.T = siz(2);       end     % # of time steps
if ~isfield(V,'Npixels'),   V.Npixels = siz(1); end     % # of pixels in ROI
if ~isfield(V,'dt'),                                    % frame duration
    fr = input('\nwhat was the frame rate for this movie (in Hz)?: ');
    V.dt = 1/fr;
end

% variables determined by the user
if ~isfield(V,'fast_poiss'),V.fast_poiss = 0;   end     % whether observations are Poisson
if ~isfield(V,'fast_nonlin'),   V.fast_nonlin   = 0; end
if V.fast_poiss && V.fast_nonlin,
    reply = input('\ncan be nonlinear observations and poisson, \ntype 1 for nonlin, 2 for poisson, anything else for neither: ');
    if reply==1,        V.fast_poiss = 0;   V.fast_nonlin = 1;
    elseif reply==2,    V.fast_poiss = 1;   V.fast_nonlin = 0;
    else                V.fast_poiss = 0;   V.fast_nonlin = 0;
    end
end
if ~isfield(V,'fast_iter_max'), V.fast_iter_max=1; end % max # of iterations before convergence

% things that matter if we are iterating to estimate parameters
if V.fast_iter_max>1;
    if V.fast_poiss || V.fast_nonlin,
        disp('\ncode does not currrently support estimating parameters for \npoisson or nonlinear observations');
        V.fast_iter_max=1;
    end

    if ~isfield(V,'fast_plot'), V.fast_plot = 0; end
    if V.fast_plot==1
        FigNum = 400;
        if V.Npixels>1, figure(FigNum), clf, end        % figure showing estimated spatial filter
        figure(FigNum+1), clf                           % figure showing estimated spike trains
        if isfield(V,'n'), siz=size(V.n); V.n(V.n==0)=NaN; if siz(1)<siz(2), V.n=V.n'; end; end
    end

    if ~isfield(V,'est_sig'),   V.est_sig   = 0; end    % whether to estimate sig
    if ~isfield(V,'est_lam'),   V.est_lam   = 0; end    % whether to estimate lam
    if ~isfield(V,'est_gam'),   V.est_gam   = 0; end    % whether to estimate gam
    if ~isfield(V,'est_a'),     V.est_a     = 0; end    % whether to estimate a
    if ~isfield(V,'est_b'),     V.est_b     = 1; end    % whether to estimate b
    if ~isfield(V,'fast_plot'), V.fast_plot = 1; end    % whether to plot results from each iteration
    if ~isfield(V,'fast_thr'),  V.fast_thr  = 0; end    % whether to threshold spike train before estimating 'a' and 'b'
    if ~isfield(V,'fast_ignore_post'), V.fast_ignore_post=0; end % whether to ignore the posterior, and just keep the last iteration
end

% normalize F if it is only a trace
if V.Npixels==1
    F=detrend(F);
    F=F-min(F)+eps;
end

%% set default model Parameters

if nargin < 3,          P       = struct;                       end
if ~isfield(P,'sig'),   P.sig   = mean(mad(F',1)*1.4826);       end
if ~isfield(P,'gam'),   P.gam   = (1-V.dt/1)*ones(V.Ncells,1);  end
if ~isfield(P,'lam'),   P.lam   = 10*ones(V.Ncells,1);          end
if ~isfield(P,'a'),     P.a     = median(F,2);                  end

if ~isfield(P,'b'),
    if V.Npixels==1, P.b = quantile(F,0.05);
    else P.b=median(F,2);
    end
end

%% define some stuff needed for est_MAP function

% for brevity and expediency
Z   = zeros(V.Ncells*V.T,1);                    % zero vector
M   = spdiags([repmat(-P.gam,V.T,1) repmat(Z,1,V.Ncells-1) (1+Z)], -V.Ncells:0,V.Ncells*V.T,V.Ncells*V.T);  % matrix transforming calcium into spikes, ie n=M*C
I   = speye(V.Ncells*V.T);                      % create out here cuz it must be reused
H1  = I;                                        % initialize memory for Hessian matrix
H2  = I;                                        % initialize memory for other part of Hessian matrix
d0  = 1:V.Ncells*V.T+1:(V.Ncells*V.T)^2;        % index of diagonal elements of TxT matrices
d1  = 1+V.Ncells:V.Ncells*V.T+1:(V.Ncells*V.T)*(V.Ncells*(V.T-1)); % index of off-diagonal elements of TxT matrices
posts = Z(1:V.fast_iter_max);                   % initialize likelihood
if numel(P.lam)==V.Ncells
    lam = V.dt*repmat(P.lam,V.T,1);             % for lik
elseif numel(P.lam)==V.Ncells*V.T
    lam = V.dt*P.lam;
else
    error('lam must either be length V.T or 1');
end

if V.fast_poiss==1
    H       = I;                                % initialize memory for Hessian matrix
    gamlnF  = gammaln(F+1);                     % for lik
    sumF    = sum(F,1)';                        % for grad & Hess
end

%% if not iterating to estimate parameters, only this is necessary
[n C posts(1)] = est_MAP(F,P);
n_best = n;
P_best = P;
V.fast_iter_tot = 1;
V.post = posts(1);
post_max = posts(1);

if V.fast_iter_max>1
    options = optimset('Display','off');        % don't show warnings for parameter estimation
    i       = 1;                                % iteration #
    i_best  = i;                                % iteration with highest likelihood
    conv    = 0;                                % whether algorithm has converged yet
else
    conv    = 1;
end

%%  if parameters are unknown, do pseudo-EM iterations
while conv == 0
    if V.fast_plot == 1, MakePlot(n,F,P,V); end % plot results from previous iteration
    i               = i+1;                      % update iteratation number
    V.fast_iter_tot = i;                        % record of total # of iterations
    P               = est_params(n,C,F,P,b);    % update parameters based on previous iteration
    [n C posts(i)]  = est_MAP(F,P);             % update inferred spike train based on new parameters

    if posts(i)>post_max || V.fast_ignore_post==1% if this is the best one, keep n and P
        n_best  = n;                            % keep n
        P_best  = P;                            % keep P
        i_best  = i;                            % keep track of which was best
        post_max= posts(i);                     % keep max posterior
    end

    % if lik doesn't change much (relatively), or returns to some previous state, stop iterating
    if  i>=V.fast_iter_max || (abs((posts(i)-posts(i-1))/posts(i))<1e-3 || any(posts(1:i-1)-posts(i))<1e-5)% abs((posts(i)-posts(i-1))/posts(i))<1e-5 || posts(i-1)-posts(i)>1e5;
        MakePlot(n,F,P,V);
        disp('convergence criteria met')
        V.post  = posts(1:i);
        conv    = 1;
    end
    sound(3*sin(linspace(0,90*pi,2000)))        % play sound to indicate iteration is over
end

V.fast_time = cputime-starttime;                % time to run code
V           = orderfields(V);                   % order fields alphabetically to they are easier to read
P_best      = orderfields(P_best);
% n_best      = n_best./repmat(max(n_best),V.T,1);

P_best.j=j;

%% fast filter function
    function [n C post] = est_MAP(F,P)

        % initialize n and C
        z = 1;                                  % weight on barrier function
        llam = reshape(1./lam',1,V.Ncells*V.T)';
        if V.fast_nonlin==1
            n = V.gauss_n;
        else
            n = 0.01+0*llam;                    % initialize spike train
        end
        C = 0*n;                                % initialize calcium
        for j=1:V.Ncells
            C(j:V.Ncells:end) = filter(1,[1, -P.gam(j)],n(j:V.Ncells:end)); %(1-P.gam(j))*P.b(j);
        end

        % precompute parameters required for evaluating and maximizing likelihood
        b           = repmat(P.b,1,V.T);       % for lik
        if V.fast_poiss==1
            suma    = sum(P.a);                 % for grad
        else
            M(d1)   = -repmat(P.gam,V.T-1,1);   % matrix transforming calcium into spikes, ie n=M*C
            ba      = P.a'*b; ba=ba(:);         % for grad
            aa      = repmat(diag(P.a'*P.a),V.T,1);% for grad
            aF      = P.a'*F; aF=aF(:);         % for grad
            e       = 1/(2*P.sig^2);            % scale of variance
            H1(d0)  = -2*e*aa;                   % for Hess
        end
        grad_lnprior  = M'*llam;                  % for grad

        % find C = argmin_{C_z} lik + prior + barrier_z
        while z>1e-13                           % this is an arbitrary threshold

            if V.fast_poiss==1
                Fexpect = P.a*(C+b')';          % expected poisson observation rate
                lik = -sum(sum(-Fexpect+ F.*log(Fexpect) - gamlnF)); % lik
            else
                if V.fast_nonlin==1
                    S = C./(C+P.k_d);
                else
                    S = C;
                end
                D = F-P.a*(reshape(S,V.Ncells,V.T))-b; % difference vector to be used in likelihood computation
                lik = e*D(:)'*D(:);             % lik
            end
            post = lik + llam'*n - z*sum(log(n));
            s    = 1;                           % step size
            d    = 1;                           % direction
            while norm(d)>5e-2 && s > 1e-3      % converge for this z (again, these thresholds are arbitrary)
                if V.fast_poiss==1
                    glik    = suma - sumF./(C+b');
                    H1(d0)  = sumF.*(C+b').^(-2); % lik contribution to Hessian
                elseif V.fast_nonlin==1
                    glik    = -2*P.a*P.k_d*D'.*(C+P.k_d).^-2;
                    H1diag  = (-P.a*P.k_d-2*(C+P.k_d).*D').*((C+P.k_d).^-4);
                    H1(d0)  = H1diag;
                else
                    glik    = -2*e*(aF-aa.*C-ba);  % gradient
                end
                g       = glik + grad_lnprior - z*M'*(n.^-1);
                H2(d0)  = n.^-2;                % log barrier part of the Hessian
                H       = H1 - z*(M'*H2*M);     % Hessian
                d   = H\g;                     % direction to step using newton-raphson
                hit = -n./(M*d);                % step within constraint boundaries
                hit=hit(hit>0);
                if any(hit<1)
                    s = min(1,0.99*min(hit));
                else
                    s = 1;
                end
                post1 = post+1;
                while post1>=post+1e-7          % make sure newton step doesn't increase objective
                    C1  = C+s*d;
                    n   = M*C1;
                    if V.fast_poiss==1
                        Fexpect = P.a*(C1+b')';
                        lik1    = -sum(sum(-Fexpect+ F.*log(Fexpect) - gamlnF));
                    else
                        if V.fast_nonlin==1
                            S1 = C1./(C1+P.k_d);
                        else
                            S1 = C1;
                        end
                        D = F-P.a*(reshape(S1,V.Ncells,V.T))-b; % difference vector to be used in likelihood computation
                        lik1 = e*D(:)'*D(:);             % lik
                    end
                    post1 = lik1 + llam'*n - z*sum(log(n));
                    s   = s/5;                  % if step increases objective function, decrease step size

                    if s<1e-20; disp('reducing s further did not increase likelihood'), break; end      % if decreasing step size just doesn't do it
                end
                C    = C1;                      % update C
                post = post1;                   % update post
            end
            z=z/10;                             % reduce z (sequence of z reductions is arbitrary)
        end

        % reshape things in the case of multiple neurons within the ROI
        n=reshape(n,V.Ncells,V.T)';
        C=reshape(C,V.Ncells,V.T)';
    end

%% Parameter Update
    function P = est_params(n,C,F,P,b)

        % generate regressor for spatial filter
        if V.est_a==1 || V.est_b==1
            if V.fast_thr==1
                CC=0*C;
                for j=1:V.Ncells
                    nsort   = sort(n(:,j));
                    nthr    = nsort(round(0.98*V.T));
                    nn      = Z(1:V.T);
                    nn(n(:,j)<=nthr)=0;
                    nn(n(:,j)>nthr)=1;
                    CC(:,j) = filter(1,[1 -P.gam(j)],nn) + (1-P.gam(j))*P.b(j);
                end
            else
                CC      = C;
            end

            if V.est_b==1
                A = [CC -1+Z(1:V.T)];
            else
                A=CC;
            end
            X = A\F';

            P.a = X(1:V.Ncells,:)';
            if V.est_b==1
                P.b = X(end,:)';
                b   = repmat(P.b,1,V.T);
            end

            D   = F-P.a*(reshape(C,V.Ncells,V.T)) - b;

            mse = D(:)'*D(:);
        end

        if V.est_a==0 && V.est_b==0 && (V.est_sig==1 || V.est_lam==1),
            D   = F-P.a*(reshape(C,V.Ncells,V.T)+b);
            mse = D(:)'*D(:);
        end

        % estimate other parameters
        if V.est_sig==1,
            P.sig = sqrt(mse)/V.T;
        end
        if V.est_lam==1,
            nnorm   = n./repmat(max(n),V.T,1);
            if numel(P.lam)==V.Ncells
                P.lam   = sum(nnorm)'/(V.T*V.dt);
                lam     = repmat(P.lam,V.T,1)*V.dt;
            else
                P.lam   = nnorm/(V.T*V.dt);
                lam     = P.lam*V.dt;
            end

        end
    end

%% MakePlot
    function MakePlot(n,F,P,V)
        if V.fast_plot == 1
            if V.Npixels>1                                     % plot spatial filter
                figure(FigNum), nrows=V.Ncells;
                for j=1:V.Ncells, subplot(1,nrows,j),
                    imagesc(reshape(P.a(:,j),V.w,V.h)),
                    title('a')
                end
            end

            figure(FigNum+1),  ncols=V.Ncells; nrows=3; END=V.T; h=zeros(V.Ncells,2);
            for j=1:V.Ncells                                  % plot inferred spike train
                h(j,1)=subplot(nrows,ncols,(j-1)*ncols+1); cla
                if V.Npixels>1, Ftemp=mean(F); else Ftemp=F; end
                plot(z1(Ftemp(2:END))+1), hold on,
                bar(z1(n_best(2:END,j)))
                title(['best iteration ' num2str(i_best)]),
                axis('tight')
                set(gca,'XTickLabel',[],'YTickLabel',[])

                h(j,2)=subplot(nrows,ncols,(j-1)*ncols+2); cla
                bar(z1(n(2:END,j)))
                if isfield(V,'n'), hold on,
                    for k=1:V.Ncells
                        stem(V.n(2:END,k)+k/10,'LineStyle','none','Marker','v','MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',2)
                    end
                end
                set(gca,'XTickLabel',[],'YTickLabel',[])
                title(['current iteration ' num2str(i)]),
                axis('tight')
            end

            subplot(nrows,ncols,j*nrows),
            plot(1:i,posts(1:i))    % plot record of likelihoods
            title(['max lik ' num2str(post_max,4), ',   lik ' num2str(posts(i),4)])
            set(gca,'XTick',2:i,'XTickLabel',2:i)
            drawnow
        end
    end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [rr,t,sqi] = Analyze_ABP_PPG_Waveforms(Waveform,Type,HRVparams,detectedQRS,subjectID)
%
%   Analyze_ABP_PPG_Waveforms(Waveform,Type,HRVparams,detectedQRS,subjectID)
%	OVERVIEW:
%       Analyze ABP or PPG waveform
%
%   INPUT:
%       Waveform    - matrix containing the a raw signal in each column
%       Type        - array containing the signal type of waveforms in each column:
%                     'APB' for ABP waveform
%                     'PPG' for PPG waveform
%       HRVparams   - struct of settings for HRV analysis
%       subjectID   - string to identify current subject
%
%   OUTPUT
%       rr          - (seconds) Vector containing RR interval
%       t           - (seconds) Time of the rr interval data
%       sqi         - (optional) Signal Quality Index; Requires a
%                     matrix with at least two columns. Column 1
%                     should be timestamps of each sqi measure, and
%                     Column 2 should be SQI on a scale from 0 to 1.
%
%   DEPENDENCIES & LIBRARIES:
%       PhysioNet Cardiovascular Signal Toolbox
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%
%   REFERENCE:
%   Vest et al. "An Open Source Benchmarked HRV Toolbox for Cardiovascular
%   Waveform and Interval Analysis" Physiological Measurement (In Press), 2018.
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Giulia Da Poian (giulia.dap@gmail.com) on Sep 6, 2017.
%       Dependent scripts written by various authors
%       (see functions for details)
%	COPYRIGHT (C) 2018
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NmbOfSigs = size(Waveform,2);
AnnotationFolder = strcat(HRVparams.writedata, filesep, 'Annotation', filesep);
if ~exist(AnnotationFolder, 'dir')
   mkdir(AnnotationFolder)
end
addpath(AnnotationFolder)

rr = [];
t = [];
sqi = [];

for i = 1:NmbOfSigs

    current_type = Type{i};
    switch current_type
        case 'PPG'
            % PPG Detection - qppg
            [PPGann] = qppg(Waveform(:,i),HRVparams.Fs);
            % PPG SQI
            [ppgsqi_numeric, ~, ppgsqi]= calculate_ppgsqi(PPGann,Waveform(:,i),HRVparams.Fs);
            % Write PPG  annotations
            write_ann(strcat(AnnotationFolder, subjectID),HRVparams,'ppg',PPGann);
            write_ann(strcat(AnnotationFolder, subjectID),HRVparams,'sqippg',PPGann(1:length(ppgsqi_numeric)),char(ppgsqi),ppgsqi_numeric);

            rr = diff(PPGann)./HRVparams.Fs;
            t = PPGann(2:end)./HRVparams.Fs;
            sqi = [PPGann'./HRVparams.Fs, ppgsqi_numeric'./100];

        case 'ABP'
            % ABP
            ABPann = run_wabp(Waveform(:,i));
            % ABP SQI
            ABPfeatures =  abpfeature(Waveform(:,i), ABPann, HRVparams.Fs);
            [BeatQ, ~] = jSQI(ABPfeatures, ABPann, Waveform(:,i));

            if ~isempty(detectedQRS)
                % Pulse Transit Time
                ptt = pulsetransit(detectedQRS, ABPann);
                % Plot BP vs PTT
                syst = ABPfeatures(:,2);
                if HRVparams.gen_figs
                    figure;
                    plot(syst,ptt(:,3)./HRVparams.Fs,'o');
                    xlabel('BP (mmHg)'); ylabel('PTT (s)');
                    title('Pulse Transit Time - BP vs PTT (ABP - QRS)')
                end
            end
            % Write ABP annotations
            write_ann(strcat(AnnotationFolder, subjectID),HRVparams,'abpm',ABPann);
            write_ann(strcat(AnnotationFolder, subjectID),HRVparams,'sqiabp',ABPann,'N',int8(BeatQ(:,1)));%write_ann(strcat(AnnotationFolder, subjectID),HRVparams,'sqiabp',BeatQ(:,1));

    end
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [ reoriented_data ] = reorient_chipper_data( symmetry, data_in)
%REORIENT_CHIPPER_ARGS Applies a transform to the chipper indexing
%
% This function transforms the actual data that a chipper function outputs
% based on the symmetry specified.
%
% A chipper function is a simplified function of a file reader.  It takes
% only three arguments: [minIndexInDimension1 maxIndexInDimension1],
% [minIndexInDimension2 maxIndexInDimension2], [subsampleInDimension1,
% subsampleInDimension2].  All other information such as the file handle
% and other auxilliary information must be contained within the chipper
% function.
%
% The SYMMETRY parameter is used to perform symmetry operations on the
% imported image data.  A combination of mirroring in the first dimension
% (m1), mirroring in the second dimension (m2) and transpose(tr) is used to
% accomplish the symmetry operations.  The 'symmetry' input parameter is
% defined as [ 1m m2 tr ], where each of the elements can be either 0 or 1
% to indicate whether to apply a transformation or not.
%
%   [ 0 0 0 ] - default setting; successive pixels on disk are interpreted
%               to fill the first dimension in a MATLAB array (vertical, if
%               viewed in imagesc/imshow without any reorienting)
%   [ 1 0 1 ] -  90 degrees CCW from [ 0 0 0 ]
%   [ 0 1 1 ] -  90 degrees  CW from [ 0 0 0 ]
%   [ 1 1 0 ] - 180 degrees     from [ 0 0 0 ]
%
%   [ 0 0 1 ] - transpose of [ 0 0 0 ];
%   [ 0 1 0 ] -  90 degrees CCW from [ 0 0 1 ]
%   [ 1 0 0 ] -  90 degrees  CW from [ 0 0 1 ]
%   [ 1 1 1 ] - 180 degrees     from [ 0 0 1 ]
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Transform data
reoriented_data=data_in;
if symmetry(1)
    reoriented_data=flipud(reoriented_data);
end
if symmetry(2)
    reoriented_data=fliplr(reoriented_data);
end
if symmetry(3)
    reoriented_data=permute(reoriented_data,[2 1 3:ndims(reoriented_data)]);
end

end

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

author: ngageoint


text:
function A = shiftAnglesFromMinus180To180(InA)

    %dont do anything to angles that are already in range
    %if we leave out this step, we get +180 shifted to -180
    idx = find(InA>180 | InA <-180);
    if(isempty(idx))
        A = InA;
    else
        A = InA;
        A(idx) = mod(InA(idx)+180,2*180)-180;  % input is in degrees
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function dump_error(list, errimg, errdet, misimg, misdet, gt, root)
% dump_error(list, error_img, err_detail, miss_image, miss_detail, 'error_dump/')
    if nargin < 7
        dump = false;
    else
        dump = true;
    end
    all_neg = errimg; %union(errimg, misimg);
    for i = 1 : length(all_neg)
        id = all_neg(i);
        fprintf('%d...\n', id);
        img = imread(list{id});
        gtt = gt{i};
        gtt = [gtt(:,[1 2]) gtt(:,3)-gtt(:,1) gtt(:,4)-gtt(:,2)];
        img = insertShape(img, 'rectangle', gtt, 'Color', 'yellow', 'LineWidth', 2);
        if(any(errimg==id))
            err = errdet{id};
            err = [err(:,[1 2]) err(:,3)-err(:,1) err(:,4)-err(:,2)];
            img = insertShape(img, 'rectangle', err, 'Color', 'yellow', 'LineWidth', 2);
        end
        if(any(misimg==id))
            err = misdet{id};
            err = [err(:,[1 2]) err(:,3)-err(:,1) err(:,4)-err(:,2)];
            img = insertShape(img, 'rectangle', err, 'Color', 'red', 'LineWidth', 2);
        end
        hold off;
        imshow(img);
        if dump
            imwrite(img, fullfile(root, [num2str(id) '.jpg']));
        end
        waitforbuttonpress;
    end

end

author: liuyuisanai

repo: RSA-for-object-detection


text:
function [err_N, sG] = ptv_local_nuclear_metric(volfix, voldef, pix_resolution, mask, singular_coefs, P_SZ, norm_type)
%% todo:
% * volfix
% * zero mean, choose std???
% * making zero mean changes gradient!

%     norm_type = 3;
    abs_type = 0;
    deps = 1e-3;

    volsz = [size(voldef, 1), size(voldef, 2), size(voldef, 3)];
    Nd = 3;
    if volsz(3) == 1
        Nd = 2;
    end

    P_SZ = round(P_SZ);
    Nimgs = size(voldef, 5);
    if Nd == 2
        P_SZ = [P_SZ, 1];
    end

    err_N = 0;
    sG = zeros([volsz, Nimgs], 'like', voldef);
    if Nd == 2
        for i1 = 1 : P_SZ(1) : volsz(1)
        for i2 = 1 : P_SZ(2) : volsz(2)
            i11 = min(i1 + P_SZ(1) - 1, volsz(1));
            i22 = min(i2 + P_SZ(2) - 1, volsz(2));
            if i11-i1 < 2 || i22 - i2 < 2
                continue;
            end
            patch = voldef(i1:i11, i2:i22, 1, :, :);
            if norm_type == 1
                patch = patch - sum(sum(sum(patch,1),2),3) / (size(patch,1)*size(patch,2)*size(patch,3));
            elseif norm_type == 2
                patch = patch - mean(patch(:));
            elseif norm_type == 3
                patch = patch - sum(sum( patch, 4), 5) / (size(patch,4)*size(patch,5));
            elseif norm_type == 4
                patch = patch - sum(sum(sum(patch,1),2),3) / (size(patch,1)*size(patch,2)*size(patch,3));
                patch = patch - sum(sum( patch, 4), 5) / (size(patch,4)*size(patch,5));
            end

            patch0 = reshape(patch, [size(patch, 1), size(patch, 2), size(patch, 3), Nimgs]);
            if abs_type == 1
                patch = abs(patch);
            elseif abs_type == 2
                patch = sqrt(patch.^2 + deps);
            end

            if ~isempty(mask)
                mask_cur = logical(mask(i1:i11, i2:i22));
            else
                mask_cur = [];
            end
            if isempty(mask_cur) || any(mask_cur(:))
                [terr, tgrad] = ptv_nuclear_metric([], patch, pix_resolution, mask_cur, singular_coefs);
                err_N = err_N + terr;

                if abs_type == 1
                    tgrad = tgrad.*sign(patch0);
                elseif abs_type == 2
                    tgrad = tgrad .* patch0 ./ sqrt(patch0.^2 + deps);
                end

                if norm_type == 1
                    tgrad = tgrad - sum(sum(sum(tgrad, 1), 2), 3) / (size(patch,1)*size(patch,2)*size(patch,3));
                elseif norm_type == 2
                    tgrad = tgrad - mean(tgrad(:));
                elseif norm_type == 3
                    tgrad = tgrad - sum(sum( tgrad, 4), 5) / (size(patch,4)*size(patch,5));
                elseif norm_type == 4
                    tgrad = tgrad - sum(sum( tgrad, 4), 5) / (size(patch,4)*size(patch,5));
                    tgrad = tgrad - sum(sum(sum(tgrad, 1), 2), 3) / (size(patch,1)*size(patch,2)*size(patch,3));
                end

                sG(i1:i11, i2:i22, 1, :) = tgrad;
            end
        end
        end
    else %3d
        for i1 = 1 : P_SZ(1) : volsz(1)
        for i2 = 1 : P_SZ(2) : volsz(2)
        for i3 = 1 : P_SZ(3) : volsz(3)
            i11 = min(i1 + P_SZ(1) - 1, volsz(1));
            i22 = min(i2 + P_SZ(2) - 1, volsz(2));
            i33 = min(i3 + P_SZ(3) - 1, volsz(3));
            if i11-i1 < 2 || i22 - i2 < 2 || i33 - i3 < 2
                continue;
            end
            patch = voldef(i1:i11, i2:i22, i3:i33, :, :);
            if norm_type == 1
                patch = patch - sum(sum(sum(patch,1),2),3) / (size(patch,1)*size(patch,2)*size(patch,3));
            elseif norm_type == 2
                patch = patch - mean(patch(:));
            elseif norm_type == 3
                patch = patch - sum(sum( patch, 4), 5) / (size(patch,4)*size(patch,5));
            elseif norm_type == 4
                patch = patch - sum(sum(sum(patch,1),2),3) / (size(patch,1)*size(patch,2)*size(patch,3));
                patch = patch - sum(sum( patch, 4), 5) / (size(patch,4)*size(patch,5));
            end

            patch0 = reshape(patch, [size(patch, 1), size(patch, 2), size(patch, 3), Nimgs]);
            if abs_type == 1
                patch = abs(patch);
            elseif abs_type == 2
                patch = sqrt(patch.^2 + deps);
            end

            if ~isempty(mask)
                mask_cur = logical(mask(i1:i11, i2:i22, i3:i33));
            else
                mask_cur = [];
            end
            if isempty(mask_cur) || any(mask_cur(:))
                [terr, tgrad] = ptv_nuclear_metric([], patch, pix_resolution, mask_cur, singular_coefs);
                err_N = err_N + terr;
                if abs_type == 1
                    tgrad = tgrad.*sign(patch0);
                elseif abs_type == 2
                    tgrad = tgrad .* patch0 ./ sqrt(patch0.^2 + deps);
                end
                if norm_type == 1
                    tgrad = tgrad - sum(sum(sum(tgrad, 1), 2), 3) / (size(patch,1)*size(patch,2)*size(patch,3));
                elseif norm_type == 2
                    tgrad = tgrad - mean(tgrad(:));
                elseif norm_type == 3
                    tgrad = tgrad - sum(sum( tgrad, 4), 5) / (size(patch,4)*size(patch,5));
                elseif norm_type == 4
                    tgrad = tgrad - sum(sum( tgrad, 4), 5) / (size(patch,4)*size(patch,5));
                    tgrad = tgrad - sum(sum(sum(tgrad, 1), 2), 3) / (size(patch,1)*size(patch,2)*size(patch,3));
                end
                sG(i1:i11, i2:i22, i3:i33, :) = tgrad;
            end
        end
        end
        end
    end
    err_N = err_N * prod(pix_resolution);
    sG = sG(1:volsz(1), 1:volsz(2), 1:volsz(3), :, :) * prod(pix_resolution);
end

% function [err_N, sG] = ptv_local_nuclear_metric(volfix, voldef, pix_resolution, mask, singular_coefs, P_SZ)
% %% todo:
% % * volfix
% % * zero mean, choose std???
% % * making zero mean changes gradient!
%
% %     P_SZ = [15,15,15];
%     volsz = [size(voldef, 1), size(voldef, 2), size(voldef, 3)];
%     Nd = 3;
%     if volsz(3) == 1
%         Nd = 2;
%     end
%
%     P_SZ = round(P_SZ);
%     Nimgs = size(voldef, 5);
%     if Nd == 2
%         P_SZ = [P_SZ, 1];
%     end
%
%     volsz_p = ceil(volsz./P_SZ) .* P_SZ;
%
%     if Nd == 2
%         volsz_p(3) = 1;
%     end
%
%     voldef_p = zeros([volsz_p, 1, Nimgs]);
%     voldef_p(1:volsz(1), 1:volsz(2), 1:volsz(3), :, :) = voldef;
%
%     if ~isempty(mask)
%         if Nd == 3
%             mask_p = zeros([volsz_p]);
%             mask_p(1:volsz(1), 1:volsz(2), 1:volsz(3)) = mask;
%         else
%             mask_p = zeros([volsz_p]);
%             mask_p(1:volsz(1), 1:volsz(2), 1) = mask;
%         end
%     end
%
%     err_N = 0;
%     sG = zeros([volsz, Nimgs], 'like', voldef);
%     if Nd == 2
%         for i1 = 1 : P_SZ(1) : volsz(1)
%         for i2 = 1 : P_SZ(2) : volsz(2)
%             patch = voldef_p(i1:i1+P_SZ(1)-1, i2:i2+P_SZ(2)-1, 1, :, :);
%             patch = patch - sum(sum(sum(patch,1),2),3) / numel(patch);
% %             patch = patch ./ std(std(std(patch,[],1),[],2),[],3);
%             if ~isempty(mask)
%                 mask_cur = logical(mask_p(i1:i1+P_SZ(1)-1, i2:i2+P_SZ(2)-1));
%             else
%                 mask_cur = [];
%             end
%             if isempty(mask_cur) || any(mask_cur(:))
%                 [terr, tgrad] = ptv_nuclear_metric([], patch, pix_resolution, mask_cur, singular_coefs);
%                 err_N = err_N + terr;
%                 tgrad = tgrad - sum(sum(sum(tgrad, 1), 2), 3) / numel(patch);
%                 sG(i1:i1+P_SZ(1)-1, i2:i2+P_SZ(2)-1, 1, :) = tgrad;
%             end
%         end
%         end
%     else %3d
%         for i1 = 1 : P_SZ(1) : volsz(1)
%         for i2 = 1 : P_SZ(2) : volsz(2)
%         for i3 = 1 : P_SZ(3) : volsz(3)
%             patch = voldef_p(i1:i1+P_SZ(1)-1, i2:i2+P_SZ(2)-1, i3:i3+P_SZ(3)-1, :, :);
%             patch = patch - sum(sum(sum(patch,1),2),3) / numel(patch);
%             if ~isempty(mask)
%                 mask_cur = logical(mask_p(i1:i1+P_SZ(1)-1, i2:i2+P_SZ(2)-1, i3:i3+P_SZ(3)-1));
%             else
%                 mask_cur = [];
%             end
%             if isempty(mask_cur) || any(mask_cur(:))
%                 [terr, tgrad] = ptv_nuclear_metric([], patch, pix_resolution, mask_cur, singular_coefs);
%                 err_N = err_N + terr;
%                 sG(i1:i1+P_SZ(1)-1, i2:i2+P_SZ(2)-1, i3:i3+P_SZ(3)-1, :) = tgrad;
%             end
%         end
%         end
%         end
%     end
%     err_N = err_N * prod(pix_resolution);
%     sG = sG(1:volsz(1), 1:volsz(2), 1:volsz(3), :, :) * prod(pix_resolution);
% end

repo: pTVreg


text:
function expInfo = exportInfo(guifile)
%EXPORTINFO     Extract useful info from a CHEBGUI object for exporting.
%
% Calling sequence
%
%   EXPINFO = EXPORTINFO(GUIFILE)
%
% where
%
%   GUIFILE:    A CHEBGUI object
%   EXPINFO:    A struct, containing fields with information for exporting to an
%               .m-file.

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

% Extract information from the GUI fields
dom = guifile.domain;
deInput = guifile.DE;
bcInput = guifile.BC;
initInput = guifile.init;

% Wrap all input strings in a cell (if they're not a cell already)
if ( isa(deInput, 'char') )
     deInput = cellstr(deInput);
end

if ( isa(bcInput, 'char') )
     bcInput = cellstr(bcInput);
end

if ( isa(initInput, 'char') )
     initInput = cellstr(initInput);
end

% Obtain useful strings describing the differential equation part:
[deString, allVarString, indVarNameDE, dummy, dummy, dummy, allVarNames] = ...
    setupFields(guifile, deInput, 'DE');

% Do some error checking before we do further printing. Check that independent
% variable name match.

% Ensure that no EIG specific variables appear:
eigNames = {'lambda'; 'lam'; 'l'};
eigMatch = zeros(size(allVarNames));

for eigNameCounter = 1:length(eigNames);
    eigMatch = eigMatch + strcmp(allVarNames, eigNames{eigNameCounter});
end

assert(~any(eigMatch), 'CHEBFUN:CHEBGUIEXPORTER:exportInfo:eig', ...
    ['Problem appears to be an eigenvalue problem. \n Please make sure ' ...
    '''l'', ''lam'' or ''lambda'' do not appear in appear in input.']);

% Obtain the independent variable name appearing in the initial condition:
useLatest = strcmpi(initInput{1}, 'Using latest solution');
if ( ~isempty(initInput{1}) && ~useLatest )
    [dummy, dummy, indVarNameInit] = ...
        setupFields(guifile, initInput, 'INIT', allVarString);
else
    indVarNameInit = {''};
end

% Make sure we don't have a discrepency in indVarNames
if ( ~isempty(indVarNameInit{1}) && ~isempty(indVarNameDE{1}) )
    if ( strcmp(indVarNameDE{1}, indVarNameInit{1}) )
        indVarNameSpace = indVarNameDE{1};
    else
        error('CHEBFUN:CHEBGUIEXPORTERBVP:exportInfo:SolveGUIbvp', 'Independent variable names do not agree')
    end
elseif ( ~isempty(indVarNameInit{1}) && isempty(indVarNameDE{1}) )
    indVarNameSpace = indVarNameInit{1};
elseif ( isempty(indVarNameInit{1}) && ~isempty(indVarNameDE{1}) )
    indVarNameSpace = indVarNameDE{1};
else
    indVarNameSpace = 't'; % Default value for IVPs
end

% Replace the 'DUMMYSPACE' variable in the DE field
deString = strrep(deString, 'DUMMYSPACE', indVarNameSpace);
deString = chebguiExporter.prettyPrintFevalString(deString, allVarNames);

% Do we want to solve the problem globally, or with timestepping?
timeSteppingSolver = ~isempty(strfind(guifile.options.ivpSolver, 'ode'));

% Add spaces to DOM and ALLVARSTRING so it looks nices once we export
dom = strrep(dom, ',', ', ');
allVarString = strrep(allVarString, ',', ', ');

%% Fill up the expInfo struct
expInfo.dom = dom;
expInfo.deInput = deInput;
expInfo.bcInput = bcInput;
expInfo.initInput = initInput;
expInfo.deString = deString;
expInfo.allVarString = allVarString;
expInfo.allVarNames = allVarNames;
expInfo.numVars = length(allVarNames);
expInfo.indVarNameSpace = indVarNameSpace;
expInfo.useLatest = useLatest;
expInfo.timeSteppingSolver = timeSteppingSolver;

%Information related to options set-up
expInfo.tol = guifile.tol;
expInfo.dampingOn = guifile.options.damping;
expInfo.discretization = guifile.options.discretization;
expInfo.plotting = guifile.options.plotting;
expInfo.ivpSolver = guifile.options.ivpSolver;
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function [f_lp, f_ms, cands, start_ths, end_ths] = full_cands_from_hiers(lps,ms,ths,n_cands)

n_hiers  = length(ms);
n_r_cand = size(n_cands,1);
assert(size(n_cands,2)==n_hiers)

% Scan all hierarchies
all_cands = cell(n_hiers,n_r_cand);
for ii=1:n_hiers
    if size(ms{ii},1)>0
        n_r_hier = ms{ii}(end,end);
        assert(length(unique(lps(:,:,ii)))==ms{ii}(1,end)-1)
        assert(ms{ii}(1,end)-1+size(ms{ii},1)==n_r_hier)

        % Get all pairs of neighboring leave regions
        [~, idx_neighbors] = seg2gridbmap(lps(:,:,ii));
        K = max(idx_neighbors.matrix_max(:)) + 1;
        neigh_pairs = unique(idx_neighbors.matrix_min+K*idx_neighbors.matrix_max);
        neigh_pairs(neigh_pairs==0) = [];
        neigh_pairs_min = mod(neigh_pairs,K);
        neigh_pairs_max = (neigh_pairs-neigh_pairs_min)/K;

        if isrow(neigh_pairs_min)
            neigh_pairs_min = neigh_pairs_min';
        end
        if isrow(neigh_pairs_max)
            neigh_pairs_max = neigh_pairs_max';
        end

        % Get the 'n_cands' top candidates from each hierarchy
        % Singletons
        if n_r_hier<=n_cands(1,ii)
            all_cands{ii,1} = (1:n_r_hier)';
        else
            all_cands{ii,1} = (n_r_hier-n_cands(1,ii)+1:n_r_hier)';
        end

    %     % Pairs, triplets
    %     [all_cands{ii,2}, all_cands{ii,3}] = ...
    %         mex_get_tree_cands(double(lps(:,:,ii))-1, double(ms{ii})-1,...
    %                              neigh_pairs_min-1, neigh_pairs_max-1,...
    %                              [n_cands(2,ii), n_cands(3,ii)]);

        % Pairs, triplets, etc.
        all_cands(ii,2:n_r_cand) = ...
        mex_get_tree_cands(double(lps(:,:,ii))-1, double(ms{ii})-1,...
                             neigh_pairs_min-1, neigh_pairs_max-1,...
                             n_cands(2:end,ii));
    else
        all_cands{1} = 1;
    end
end

% Put all them together (for each hierarchy)
full_cands = cell(n_hiers,1);
for ii=1:n_hiers
    % Pre-allocate
    n_tot_cand = 0;
    for jj=1:n_r_cand
        n_tot_cand = n_tot_cand + size(all_cands{ii,jj},1);
    end
    full_cands{ii} = zeros(n_tot_cand,n_r_cand);

    if n_tot_cand
        % Store
        curr_n = 0;
        for jj=1:n_r_cand
            if size(all_cands{ii,jj},1)>0
                full_cands{ii}(curr_n+1:curr_n+size(all_cands{ii,jj},1),1:jj) = all_cands{ii,jj};
                curr_n = curr_n+size(all_cands{ii,jj},1);
            end
        end
    end
end

% Compute the 'unique' hierarchy
lps = double(lps);
luts = cell(n_hiers,1);
for ii=1:n_hiers
    needed_regs = unique(full_cands{ii}(:));
    needed_regs(needed_regs==0) = [];

    if isempty(needed_regs) || isempty(ms{ii})
        lps(:,:,ii) = ones(size(lps,1),size(lps,2));
    else
        [lps(:,:,ii),ms{ii},luts{ii}] = mex_prune_tree_to_regions(lps(:,:,ii)-1,ms{ii}-1,needed_regs-1);
    end
end

% Fuse
[f_lp, f_ms] = fuse_bpts(lps,ms);

% Reduce zeros
f_ms = f_ms(:,logical(sum(f_ms)>0));

% Redo cands (relabel)
n_int = length(unique(f_lp));
curr_n_regs = n_int;
cands = [];
for jj=1:n_hiers
    if ~isempty(ms{jj})
        curr_part = lps(:,:,jj);
        assert(max(curr_part(:))+size(ms{jj},1)==size(luts{jj},1))
        new_cands = zeros(size(full_cands{jj}));
        for xx=1:size(full_cands{jj},1)
            for yy=1:size(full_cands{jj},2)
                if full_cands{jj}(xx,yy)>0
                    pos = logical(full_cands{jj}(xx,yy)==luts{jj}(:,2));
                    if sum(pos)~=1
                        error('Oh oh')
                    end
                    new_cands(xx,yy) = curr_n_regs+luts{jj}(pos,1);
                end
            end
        end
        cands = [cands; new_cands]; %#ok<AGROW>
        curr_n_regs = curr_n_regs + ms{jj}(end,end);
        % Sanity check
        assert(ms{jj}(end,end)==length(unique(lps(:,:,jj))) + size(ms{jj},1))
    else
        if ~isempty(full_cands{jj})
            cands = [cands; full_cands{jj}]; %#ok<AGROW>
        end
    end
end

% Sanity checks
if ~isempty(cands)
    assert(isequal(size(cell2mat(full_cands)), size(cands)))
else
    assert(isempty(cell2mat(full_cands)))
end

% Redo ths
start_ths = zeros(1,n_int);
end_ths = zeros(1,n_int);
for jj=1:n_hiers
    if ~isempty(ms{jj})
        curr_part = lps(:,:,jj);
        assert(max(curr_part(:))+size(ms{jj},1)==size(luts{jj},1))

        start_ths = [start_ths ths{jj}.start_ths(luts{jj}(:,2))]; %#ok<AGROW>
        end_ths   = [end_ths ths{jj}.end_ths(luts{jj}(:,2))]; %#ok<AGROW>
    end
end
assert(length(start_ths)==n_int+size(f_ms,1))

end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function [varargout] = istable(varargin)

% ISTABLE is a drop-in replacement for the same function that was
% introduced in MATLAB R2013b.
%
% In all MATLAB versions prior to 2013b this function returns false.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% try to automatically remove incorrect compat folders from the path, see https://github.com/fieldtrip/fieldtrip/issues/899

if isempty(strfind(mfilename('fullpath'), matlabroot))
  % automatic cleanup does not work when FieldTrip is installed inside the MATLAB path, see https://github.com/fieldtrip/fieldtrip/issues/1527

  alternatives = which(mfilename, '-all');
  if ~iscell(alternatives)
    % this is needed for octave, see https://github.com/fieldtrip/fieldtrip/pull/1171
    alternatives = {alternatives};
  end

  keep = true(size(alternatives));
  for i=1:numel(alternatives)
    keep(i) = keep(i) && ~any(alternatives{i}=='@');  % exclude methods from classes
    keep(i) = keep(i) && alternatives{i}(end)~='p';   % exclude precompiled files
  end
  alternatives = alternatives(keep);

  if exist(mfilename, 'builtin') || any(strncmp(alternatives, matlabroot, length(matlabroot)) & cellfun(@isempty, strfind(alternatives, fullfile('private', mfilename))))
    % remove this directory from the path
    p = fileparts(mfilename('fullpath'));
    warning('removing "%s" from your path, see http://bit.ly/2SPPjUS', p);
    rmpath(p);
    % call the original MATLAB function
    if exist(mfilename, 'builtin')
      [varargout{1:nargout}] = builtin(mfilename, varargin{:});
    else
      [varargout{1:nargout}] = feval(mfilename, varargin{:});
    end
    return
  end

end % automatic cleanup of compat directories

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is where the actual replacement code starts
% function tf = istable(input)

tf = false;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output arguments

varargout = {tf};

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function gIX = KmeansSubdivide(numK2,gIX,M_0)
gIX_old = gIX;
U = unique(gIX);
        for i = 1:length(U),
            IX = find(gIX_old == U(i));
            M_sub = M_0(IX,:);

            if numK2<length(IX),
                [gIX_sub,C] = kmeans(M_sub,numK2,'distance','correlation');
            else
                [gIX_sub,C] = kmeans(M_sub,length(IX),'distance','correlation');
            end
            gIX(IX) = (i-1)*numK2+gIX_sub;
        end
end

author: xiuyechen

repo: FishExplorer


%==============================================================================
% This code is part of the VAMPIRE app for the Matlab-based toolbox
%  FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR/VAMPIRE.m
%==============================================================================
%
% 2D Example for Multilevel Mass-Preserving Image Registration using VAMPIRE
%
% 	- data                 synthetic 2D Gaussian blobs (level 3:8, full
%                          resolution: 256x256)
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             SSD
%   - pre-registration     none
%   - regularizer          mfHyperElastic
%   - optimizer            Gauss-Newton with ArmijoBacktrack linesearch

setup2DGaussianData;

% prepare the plot
FAIRplots('clear')
Dshow = @(T,R,omega,m) viewImage2D(128+(T-R)/2,omega,m,'colormap',gray(256));
FAIRplots('set','Dshow',Dshow);

% initialize the regularizer for the non-parametric part
alpha       = 100;
alphaLength = 1;
alphaArea   = 0;
alphaVolume = 1;
[reg,regOptn] = regularizer('reset', 'regularizer', 'mfHyperElastic', ...
    'alpha',alpha, 'alphaLength', alphaLength, 'alphaArea', alphaArea, ...
    'alphaVolume', alphaVolume);

% finally: run the Mass-Preserving Non-Parametric Image Registration
NPIRpara            = optPara('NPIR-GN');
NPIRpara.lineSearch = @ArmijoDiffeomorphic;
NPIRpara.solver     = @VAMPIREsolveGN_PCG;

[yc,wc,his] = MLIR(ML, 'NPIRobj', @VAMPIRENPIRobjFctn, ...
    'parametric', false, 'NPIRpara', NPIRpara, 'minLevel', 5);

% [reg,regOptn] = regularizer('reset','regularizer','mfElastic',...
%   'alpha',alpha,'alphaLength',alphaLength,'alphaArea',alphaArea,...
%   'alphaVolume',alphaVolume);
%
% % finally: run the  Non-Parametric Image Registration
% [yc,wc,his] = MLIR(ML, 'NPIRobj', @NPIRobjFctn, ...
%   'parametric', false, 'minLevel', 4);

%% Plot Results
% Compute resulting image: dataT(yc) * det(D(yc))
Topt = reshape(linearInter(dataT,omega,center(yc,m)) .* geometry(yc,m,'Jac','omega',omega), m);
figure;
subplot(2,2,1);
viewImage2D(dataT,omega,m,'colormap', 'gray(256)');
hold on; plotGrid(center(yc, m), omega, m, 'spacing', [5 5]); axis off; hold off;
title('Template (T) & Grid')
subplot(2,2,2);
viewImage2D(dataR,omega,m,'colormap', 'gray(256)'); axis off;
title('Reference (R)')
subplot(2,2,3);
viewImage2D(Topt,omega,m,'colormap', 'gray(256)'); axis off;
title('VAMPIRE result (Topt)')
subplot(2,2,4);
viewImage2D(abs(Topt-dataR),omega,m,'colormap', 'gray(256)'); axis off;
title('Absolute difference image of R and Topt');


text:
% Copyright 2018 Marc René Schädler
%
% This file is part of the mobile hearing aid prototype project
% The the mobile hearing aid prototype project 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.
%
% The mobile hearing aid prototype project 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 the mobile hearing aid prototype project. If not, see http://www.gnu.org/licenses/.

function thresholds = measure_thresholds(freqs, ear)
  thresholds = nan(size(freqs));
  for i=1:length(freqs)
    threshold = measure_sweep(freqs(i), ear);
    if ~isempty(threshold)
      thresholds(i) = threshold;
    end
  end
end

author: m-r-s

repo: hearingaid-prototype


text:
%DEMO_REGRESSION1  Regression problem demonstration for 2-input
%                  function with Gaussian process
%
%  Description
%    The regression problem consist of a data with two input
%    variables and one output variable with Gaussian noise. The
%    model constructed is following:
%
%    The observations y are assumed to satisfy
%
%         y = f + e,    where e ~ N(0, s^2)
%
%    where f is an underlying function, which we are interested in.
%    We place a zero mean Gaussian process prior for f, which
%    implies that at the observed input locations latent values
%    have prior
%
%         f ~ N(0, K),
%
%    where K is the covariance matrix, whose elements are given as
%    K_ij = k(x_i, x_j | th). The function k(x_i, x_j | th) is
%    covariance function and th its parameters.
%
%    Since both likelihood and prior are Gaussian, we obtain a
%    Gaussian marginal likelihood
%
%        p(y|th) = N(0, K + I*s^2).
%
%   By placing a prior for parameters, p(th), we can find
%   the maximum a posterior (MAP) estimate for them by maximizing
%
%       argmax   log p(y|th) + log p(th).
%         th
%
%   An approximation for the posterior of the parameters, can be
%   found using Markov chain Monte Carlo (MCMC) methods. We can
%   integrate over the parameters also with other integration
%   approximations such as grid integration.
%
%   After finding MAP estimate or posterior samples of
%   parameters, we can use them to make predictions for f_new:
%
%       p(f_new | y, th) = N(m, S),
%
%          m = K_nt*(K + I*s^2)^(-1)*y
%          S = K_new - K_nt*(K + I*s^2)^(-1)*K_tn
%
%   where K_new is the covariance matrix of new f, and K_nt between
%   new f and training f.
%
%   For more detailed discussion of Gaussian process regression see,
%   for example, Rasmussen and Williams (2006) or Vanhatalo and
%   Vehtari (2008)
%
%   The demo is organised in three parts:
%     1) data analysis with MAP estimate for the parameters
%     2) data analysis with grid integration over the parameters
%     3) data analysis with MCMC integration over the parameters
%
%  See also DEMO_*
%
%  References:
%    Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian
%    Processes for Machine Learning. The MIT Press.
%
%    Vanhatalo, J. and Vehtari, A. (2008). Modelling local and global
%    phenomena with sparse Gaussian processes. Proceedings of the 24th
%    Conference on Uncertainty in Artificial Intelligence,

% Copyright (c) 2008-2010 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari

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

%========================================================
% PART 1 data analysis with full GP model
%========================================================
disp('GP with Gaussian noise model')

% Load the data
S = which('demo_regression1');
L = strrep(S,'demo_regression1.m','demodata/dat.1');
data=load(L);
x = [data(:,1) data(:,2)];
y = data(:,3);
[n, nin] = size(x);

% Now 'x' consist of the inputs and 'y' of the output.
% 'n' and 'nin' are the number of data points and the
% dimensionality of 'x' (the number of inputs).

% ---------------------------
% --- Construct the model ---
%
% First create structures for Gaussian likelihood and squared
% exponential covariance function with ARD
lik = lik_gaussian('sigma2', 0.2^2);
gpcf = gpcf_sexp('lengthScale', [1.1 1.2], 'magnSigma2', 0.2^2)

% Set some priors
pn = prior_logunif();
lik = lik_gaussian(lik,'sigma2_prior', pn);
pl = prior_unif();
pm = prior_sqrtunif();
gpcf = gpcf_sexp(gpcf, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);

% Following lines do the same since the default type is FULL
%gp = gp_set('type','FULL','lik',lik,'cf',gpcf);
gp = gp_set('lik', lik, 'cf', gpcf);

% Demostrate how to evaluate covariance matrices.
% K contains the covariance matrix without noise variance
%  at the diagonal (the prior covariance)
% C contains the covariance matrix with noise variance at
% the diagonal (the posterior covariance)
example_x = [-1 -1 ; 0 0 ; 1 1];
[K, C] = gp_trcov(gp, example_x)

% What has happend this far is the following
% - we created structures 'gpcf' and 'lik', which describe
%   the properties of the covariance function and Gaussian likelihood (see
%   gpcf_sexp and lik_gaussian for more details)
% - we created structures that describe the prior of the length-scale
%   and magnitude of the squared exponential covariance function and
%   the prior of the noise variance. These structures were set into
%   'gpcf' and 'lik' (see prior_* for more details)
% - we created a GP structure 'gp', which has among others 'gpcf'
%   and 'lik' structures.  (see gp_set for more details)

% -----------------------------
% --- Conduct the inference ---
%
% We will make the inference first by finding a maximum a posterior
% estimate for the parameters via gradient based optimization.
% After this we will use grid integration and Markov chain Monte
% Carlo sampling to integrate over the parameters.

% --- MAP estimate  ---
disp(' MAP estimate for the parameters')
% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);

% get optimized parameter values for display
[w,s]=gp_pak(gp);
% display exp(w) and labels
disp(s), disp(exp(w))

% For last, make predictions of the underlying function on a dense
% grid and plot it. Below Eft_map is the predictive mean and
% Varf_map the predictive variance.
[xt1,xt2]=meshgrid(-1.8:0.1:1.8,-1.8:0.1:1.8);
xt=[xt1(:) xt2(:)];
[Eft_map, Varft_map] = gp_pred(gp, x, y, xt);

% Plot the prediction and data
figure(1)
clf
mesh(xt1, xt2, reshape(Eft_map,37,37));
hold on
plot3(x(:,1), x(:,2), y, '*')
axis on;
title('The predicted underlying function and the data points (MAP solution)');

% --- Grid integration ---
disp(' Grid integration over the parameters')
% Perform the grid integration and make predictions
% We can also get predictions for marginals of f, which are not
% Gaussian due to integration over the hyperparameters
[gp_array, P_TH, th, Eft_ia, Varft_ia, fx_ia, x_ia] = ...
    gp_ia(gp, x, y, xt, 'int_method', 'grid');

% Plot the predictions for two input locations
figure(2)
clf
subplot(2,1,1)
plot(x_ia(100,:), fx_ia(100,:))
title('p(f|D) at input location (-1.6, 0.7)');
subplot(2,1,2)
plot(x_ia(400,:), fx_ia(400,:))
title('p(f|D) at input location (-0.8, 1.1)');

% --- MCMC ---
disp(' MCMC integration over the parameters')
[gp_rec,g,opt] = gp_mc(gp, x, y, 'nsamples', 220,'display',20);

% After sampling we delete the burn-in and thin the sample chain
gp_rec = thin(gp_rec, 21, 2);

% Make the predictions
[Eft_mc, Varft_mc] = gp_pred(gp_rec, x, y, xt);

figure(1)
clf
subplot(1,2,1)
mesh(xt1, xt2, reshape(Eft_map,37,37));
hold on
plot3(x(:,1), x(:,2), y, '*')
axis on;
title(['The predicted underlying function ';
       'and the data points (MAP solution)']);
subplot(1,2,2)
mesh(xt1, xt2, reshape(Eft_mc,37,37));
hold on
plot3(x(:,1), x(:,2), y, '*')
axis on;
title(['The predicted underlying function  ';
       'and the data points (MCMC solution)']);
set(gcf,'pos',[93 511 1098 420])

% We can compare the posterior samples of the parameters to the MAP
% estimate that we got from optimization
figure(3)
clf
subplot(1,2,1)
plot(gp_rec.cf{1}.lengthScale)
title('The sample chain of length-scales')
subplot(1,2,2)
plot(gp_rec.cf{1}.magnSigma2)
title('The sample chain of magnitude')
set(gcf,'pos',[93 511 1098 420])

figure(4)
clf
subplot(1,4,1)
hist(gp_rec.cf{1}.lengthScale(:,1))
hold on
plot(gp.cf{1}.lengthScale(1), 0, 'rx', 'MarkerSize', 11, 'LineWidth', 2)
title('Length-scale 1')
subplot(1,4,2)
hist(gp_rec.cf{1}.lengthScale(:,2))
hold on
plot(gp.cf{1}.lengthScale(2), 0, 'rx', 'MarkerSize', 11, 'LineWidth', 2)
title('Length-scale 2')
subplot(1,4,3)
hist(gp_rec.cf{1}.magnSigma2)
hold on
plot(gp.cf{1}.magnSigma2, 0, 'rx', 'MarkerSize', 11, 'LineWidth', 2)
title('magnitude')
subplot(1,4,4)
hist(gp_rec.lik.sigma2)
hold on
plot(gp.lik.sigma2, 0, 'rx', 'MarkerSize', 11, 'LineWidth', 2)
title('Noise variance')
legend('MCMC samples', 'MAP estimate')
set(gcf,'pos',[93 511 1098 420])

% Sample from two posterior marginals and plot them alongside
% with the MAP and grid integration results
% gpmc_preds returns the predictive mean of the latent function
% with every sampled parameter value.
[Eft_mcs, Varft_mcs] = gpmc_preds(gp_rec, x, y, xt);
sf = normrnd(Eft_mcs(100,:), sqrt(Varft_mcs(100,:)));
sf2 = normrnd(Eft_mcs(400,:), sqrt(Varft_mcs(400,:)));

figure(2)
subplot(1,2,1)
[N,X] = hist(sf);
hist(sf)
hold on
plot(x_ia(100,:), max(N)/max(fx_ia(100,:))*fx_ia(100,:), 'k')
ff = norm_pdf(x_ia(100,:)', Eft_map(100), sqrt(Varft_map(100)));
plot(x_ia(100,:), max(N)/max(ff)*ff, 'r', 'lineWidth', 2)
set(gca, 'Ytick', [])
title('p(f|D) at input location (-1.6, 0.7)');
xlim([0 1])

subplot(1,2,2)
[N,X] = hist(sf2);
hist(sf2)
hold on
plot(x_ia(400,:), max(N)/max(fx_ia(400,:))*fx_ia(400,:), 'k')
ff = norm_pdf(x_ia(400,:)', Eft_map(400), sqrt(Varft_map(400)));
plot(x_ia(400,:), max(N)/max(ff)*ff, 'r', 'lineWidth', 2)
set(gca, 'Ytick', [])
title('p(f|D) at input location (-0.8, 1.1)');
xlim([-1.2 -0.6])

disp('Done')

% ========================
% Print figures for manual
% ========================
% $$$ sf = normrnd(Eft_mc(100,:), sqrt(Varft_mc(100,:)));
% $$$ sf2 = normrnd(Eft_mc(400,:), sqrt(Varft_mc(400,:)));
% $$$
% $$$ figure
% $$$ subplot(1,2,1)
% $$$ [N,X] = hist(sf);
% $$$ hist(sf)
% $$$ h = findobj(gca,'Type','patch');
% $$$ set(h,'FaceColor','w','EdgeColor','k')
% $$$ hold on
% $$$ plot(x_ia(100,:), max(N)/max(fx_ia(100,:))*fx_ia(100,:), 'k')
% $$$ ff = norm_pdf(x_ia(100,:)', Eft_map(100), sqrt(Varft_map(100)));
% $$$ plot(x_ia(100,:), max(N)/max(ff)*ff, 'k', 'lineWidth', 2)
% $$$ set(gca, 'Ytick', [])
% $$$ xlim([0 1])
% $$$ ylim([0 110])
% $$$
% $$$ subplot(1,2,2)
% $$$ [N,X] = hist(sf2);
% $$$ hist(sf2)
% $$$ h = findobj(gca,'Type','patch');
% $$$ set(h,'FaceColor','w','EdgeColor','k')
% $$$ hold on
% $$$ plot(x_ia(400,:), max(N)/max(fx_ia(400,:))*fx_ia(400,:), 'k')
% $$$ ff = norm_pdf(x_ia(400,:)', Eft_map(400), sqrt(Varft_map(400)));
% $$$ plot(x_ia(400,:), max(N)/max(ff)*ff, 'k', 'lineWidth', 2)
% $$$ set(gca, 'Ytick', [])
% $$$ xlim([-1.2 -0.5])
% $$$
% $$$
% $$$ set(gcf,'units','centimeters');
% $$$ set(gcf,'pos',[15 14 7 5])
% $$$ set(gcf,'paperunits',get(gcf,'units'))
% $$$ set(gcf,'paperpos',get(gcf,'pos'))
% $$$ print -depsc2 /proj/bayes/jpvanhat/software/doc/GPstuffDoc/pics/demo_regression1_fig3.eps
% $$$
% $$$
% $$$ figure(4)
% $$$ clf, subplot(1,4,1)
% $$$ hist(gp_rec.cf{1}.lengthScale(:,1))
% $$$ h = findobj(gca,'Type','patch');
% $$$ set(h,'FaceColor','w','EdgeColor','k')
% $$$ hold on
% $$$ plot(gp.cf{1}.lengthScale(1), 0, 'kx', 'MarkerSize', 11, 'LineWidth', 2)
% $$$ xlabel('Length-s 1')
% $$$ xlim([0.3 1.6])
% $$$
% $$$ subplot(1,4,2)
% $$$ hist(gp_rec.cf{1}.lengthScale(:,2))
% $$$ h = findobj(gca,'Type','patch');
% $$$ set(h,'FaceColor','w','EdgeColor','k')
% $$$ hold on
% $$$ plot(gp.cf{1}.lengthScale(2), 0, 'kx', 'MarkerSize', 11, 'LineWidth', 2)
% $$$ xlabel('Length-s 2')
% $$$ xlim([0.4 1.4])
% $$$
% $$$ subplot(1,4,3)
% $$$ hist(gp_rec.cf{1}.magnSigma2)
% $$$ h = findobj(gca,'Type','patch');
% $$$ set(h,'FaceColor','w','EdgeColor','k')
% $$$ hold on
% $$$ plot(gp.cf{1}.magnSigma2, 0, 'kx', 'MarkerSize', 11, 'LineWidth', 2)
% $$$ xlabel('magnitude')
% $$$ xlim([0.5 6])
% $$$
% $$$ subplot(1,4,4)
% $$$ hist(gp_rec.lik.sigma2)
% $$$ h = findobj(gca,'Type','patch');
% $$$ set(h,'FaceColor','w','EdgeColor','k')
% $$$ hold on
% $$$ plot(gp.lik.sigma2, 0, 'kx', 'MarkerSize', 11, 'LineWidth', 2)
% $$$ xlabel('Noise variance')
% $$$ xlim([0.03 0.06])
% $$$ set(gca, 'Xtick', [0.03 0.06])
% $$$
% $$$ set(gcf,'units','centimeters');
% $$$ set(gcf,'pos',[15 14 11 5])
% $$$ set(gcf,'paperunits',get(gcf,'units'))
% $$$ set(gcf,'paperpos',get(gcf,'pos'))
% $$$
% $$$ print -depsc2 /proj/bayes/jpvanhat/software/doc/GPstuffDoc/pics/demo_regression1_fig2.eps


text:
function data = fixtrialdef(data)

% FIXTRIALDEF adds a trl matrix to the raw data configuration
% which is constructed on the fly, assuming that the trials are
% consecutive segments of a continuous recording

% Copyright (C) 2009-2010, Robert Oostenveld and Jan-Mathijs Schoffelen

if isfield(data, 'sampleinfo')
  return;
end

if ~isfield(data, 'cfg')
  % fieldtrip raw data structures are expected to have a cfg
  data.cfg = [];
end

hastrial   = isfield(data, 'trial');
hastime    = isfield(data, 'time');
hasfsample = isfield(data, 'fsample');

if ~hasfsample && isfield(data, 'time')
  data.fsample = median(1./diff(data.time{1}));
end

if hastrial,
  ntrial = length(data.trial);
else
  ntrial = dimlength(data, 'rpt');
  if ~isfinite(ntrial) && strcmp(data.dimord(1:6), 'rpttap') && isfield(data, 'cumtapcnt'),
    ntrial = numel(data.cumtapcnt);
  elseif ~isfinite(ntrial)
    ntrial = 1;
  end
end

trl = ft_findcfg(data.cfg, 'trl');

nsmp = zeros(ntrial,1);
if hastrial,
  for i=1:ntrial
    nsmp(i) = size(data.trial{i}, 2);
  end
elseif ~isempty(trl)
  nsmp = trl(:,2) - trl(:,1) + 1;
end

if size(trl,1)~=numel(nsmp)
  warning_once('the trial definition in the configuration is inconsistent with the actual data');
  trl = [];
elseif isempty(trl) || ~all(nsmp==trl(:,2)-trl(:,1)+1)
  warning_once('the data does not contain a trial definition, assuming that the trials are consecutive segments of a continuous recording');
  % construct a trial definition on the fly, assume that the trials are
  % consecutive segments of a continuous recording
  if ntrial==1,
    begsample = 1;
  else
    begsample = cat(1, 0, cumsum(nsmp(1:end-1))) + 1;
  end
  endsample = begsample + nsmp - 1;

  offset    = zeros(ntrial,1);
  if hastime,
    for i=1:ntrial
      offset(i) = time2offset(data.time{i}, data.fsample);
    end
  end
  trl = [begsample endsample offset];

elseif size(trl,1)~=ntrial
  warning_once('the trial definition in the configuration is inconsistent with the actual data');
  trl = [];
elseif nsmp~=(trl(:,2)-trl(:,1)+1)
  warning_once('the trial definition in the configuration is inconsistent with the actual data');
  trl = [];
end

if ~isfield(data, 'sampleinfo') && ~isempty(trl)
  data.sampleinfo = trl(:, 1:2);
elseif ~isfield(data, 'sampleinfo') && isempty(trl)
  warning_once('failed to create sampleinfo field');
end

if (~isfield(data, 'trialinfo') || isempty(data.trialinfo)) && ~isempty(trl) && size(trl, 2) > 3,
  data.trialinfo = trl(:, 4:end);
end

% if data is not raw then it does not make sense to keep the sampleinfo
if ~hastrial && isfield(data, 'sampleinfo')
  data = rmfield(data, 'sampleinfo');
end

author: PatternRecognition

repo: OpenBMI


text:
%  lyapunov_equation solves the equation V=T*V*T'+Q
%
%  ::
%
%    [V,retcode]=lyapunov_equation(T,Q)
%    [V,retcode]=lyapunov_equation(T,Q,options)
%
%  Args:
%     - T :
%     - Q :
%     - options :
%
%  Returns:
%     :
%     - V :
%     - retcode :
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
close all;
clearvars;
clc;
rng default;

% Ambient space dimension
M = 10;
% Number of subspaces
K = 4;
% common dimension for each subspace
D = 5;

% Construct bases for random subspaces
bases = spx.data.synthetic.subspaces.random_subspaces(M, K, D);
cluster_sizes = [10 4 4 8];
data_points = spx.data.synthetic.subspaces.uniform_points_on_subspaces(bases, cluster_sizes);

author: indigits

repo: sparse-plex


text:
function varargout = gammaincinv_x(varargin)

switch class(varargin{1})

    case 'double'
        varargout{1} = gammaincinv(varargin{1},varargin{2});

    case 'char'

        operator = CreateBasicOperator('increasing','positive','callback');
        operator.range = [0 inf];
        operator.domain = [0 1];

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

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

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [StencilC, CenterC] = stencilCrop(Stencil, Center)
%------------------------------------------------------------------------------
%
% This function crops a stencil at its bounds, in case all values are
% zero on the bounds and the center remains within the stencil.
% For dimensions <=3 no (more) cropping is performed.
%
% Stencil  = input stencil
%
% Center   = center of Stencil
%
% StencilC = output stencil (cropped)
%
% CenterC  = center of StencilC
%
% See also: stencilR2Q
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: January 24, 2001.
% (c) 1999-2002 Stichting CWI, Amsterdam.
%------------------------------------------------------------------------------
o=[0 0];
if ~all(size(o) == size(Center))
  error(' stencilCrop - unexpected dimensions of Center ')
else
  clear o;
end
StencilC = Stencil;
CenterC = Center;
%
[n, m] = size(Stencil);
for j = n:-1:4
   [nC, mC] = size(StencilC);
   v = StencilC(nC,:);
   if any(v) || CenterC(1)>=nC
     break;
   else
     StencilC = StencilC(1:(nC-1),:);
   end
end
for j = m:-1:4
   [nC, mC] = size(StencilC);
   v = StencilC(:,mC)';
   if any(v) || CenterC(2)>=mC
     break;
   else
     StencilC = StencilC(:,1:(mC-1));
   end
end
%
[n, m] = size(StencilC);
for j = 1:(n-3)
   [nC, mC] = size(StencilC);
   v = StencilC(1,:);
   if any(v) || CenterC(1)<=1
     break;
   else
     StencilC = StencilC(2:nC,:);
     CenterC = [(CenterC(1)-1) CenterC(2)];
   end
end
for j = 1:(m-3)
   [nC, mC] = size(StencilC);
   v = StencilC(:,1)';
   if any(v) || CenterC(2)<=1
     break;
   else
     StencilC = StencilC(:,2:mC);
     CenterC = [CenterC(1) (CenterC(2)-1)];
   end
end
%


text:
function varargout = spm_bms_display(BMS,action)
% Display results from BMS Maps
% FORMAT spm_bms_display(BMS,action)
%
% Input:
% BMS    - BMS containing details of excursion set
% action - 'Init' (Initialise)
%          'do_plot' (plot voxel results)
%          'save' (save results as NIfTI image)
%          'overlays' (options overlays menu)
%__________________________________________________________________________
% Copyright (C) 2009-2019 Wellcome Trust Centre for Neuroimaging

% Maria Joao Rosa
% $Id: spm_bms_display.m 7577 2019-04-24 08:59:56Z guillaume $

% Main options (action)
% =========================================================================
switch action

    % Inititalise - action: 'Init'
    % =====================================================================
    case 'Init'

    % Initialise variables
    % ---------------------------------------------------------------------
    xSPM  = BMS.xSPM;
    M     = xSPM.xVol.M;
    iM    = xSPM.iM;
    DIM   = xSPM.xVol.DIM;

    try
        if strcmp(spm('CheckModality'),'EEG')
            datatype = {...
                'Volumetric (2D/3D)',...
                'Scalp-Time',...
                'Scalp-Frequency',...
                'Time-Frequency',...
                'Frequency-Frequency'};
            selected = spm_input('Data Type: ','+1','m',datatype);
            datatype = datatype{selected};
        else
            datatype = 'Volumetric (2D/3D)';
        end
    catch
        datatype     = 'Volumetric (2D/3D)';
    end

    switch datatype
        case 'Volumetric (2D/3D)'
            units    = {'mm' 'mm' 'mm'};
        case 'Scalp-Time'
            units    = {'mm' 'mm' 'ms'};
        case 'Scalp-Frequency'
            units    = {'mm' 'mm' 'Hz'};
        case 'Time-Frequency'
            units    = {'Hz' 'ms' ''};
        case 'Frequency-Frequency'
            units    = {'Hz' 'Hz' ''};
        otherwise
            error('Unknown data type.');
    end

    title = 'Bayesian Model Selection';
    str   = xSPM.str;

    % Initialise figures
    % ---------------------------------------------------------------------
    [Finter,Fgraph,CmdLine] = spm('FnUIsetup','BMS: Results');
    FS = spm('FontSizes');
    WS = spm('WinScale');
    PF = spm_platform('fonts');

    % Clear satellite figure if it exists
    % ---------------------------------------------------------------------
    hSat = findobj('tag','Satellite');
    spm_figure('clear',hSat);

    % Setup Finter (interactive window)
    % ---------------------------------------------------------------------
    spm_figure('Clear',Finter);
    spm('FigName','BMS results',Finter,CmdLine);
    Finter  = spm_figure('GetWin',Finter);

    hReg    = uicontrol(Finter,'Style','Frame','Position',...
        [001 001 400 190].*WS,'BackgroundColor',spm('Colour'));

    hFResUi = uicontrol(Finter,'Style','Frame','Position',...
        [008 007 387 178].*WS);

    [hReg,xyz] = spm_XYZreg('InitReg',hReg,M,DIM,[0;0;0]);

    % Draw MIP
    % ---------------------------------------------------------------------
    hMIPax = axes('Parent',Fgraph,'Position',...
        [0.125 0.5450 0.59 0.40],'Visible','off');
    hMIPax = spm_mip_ui(xSPM.Z,xSPM.XYZmm,M,DIM,hMIPax,units);
    spm_XYZreg('XReg',hReg,hMIPax,'spm_mip_ui');
    hTitAx = axes('Parent',Fgraph,'Position',[0.02 0.95 0.86 0.02],...
        'Visible','off');
    text(0.5,0,title,'Parent',hTitAx,'HorizontalAlignment','center',...
        'VerticalAlignment','baseline','FontWeight','Bold','FontSize',FS(14))

    text(240,260,str,'Interpreter','TeX','FontSize',FS(14),...
        'Fontweight','Bold','Parent',hMIPax)

    % Print BMSresults: Results directory & thresholding info
    %----------------------------------------------------------------------
    hResAx = axes('Parent',Fgraph,...
        'Position',[0.160 0.510 0.45 0.05],...
        'DefaultTextVerticalAlignment','baseline',...
        'DefaultTextFontSize',FS(9),...
        'DefaultTextColor',[1,1,1]*.7,...
        'Units','points',...
        'Visible','off');
    AxPos = get(hResAx,'Position'); set(hResAx,'YLim',[0,AxPos(4)])
    h     = text(0,24,'BMSresults:','Parent',hResAx,...
        'FontWeight','Bold','FontSize',FS(14));
    text(get(h,'Extent')*[0;0;1;0],24,spm_file(pwd,'short30'),'Parent',hResAx)
    text(0,12,sprintf('Threshold: %0.2d',BMS.xSPM.thres),'Parent',hResAx)

    % Store handles of results section Graphics window objects
    %----------------------------------------------------------------------
    H  = get(Fgraph,'Children');
    H  = findobj(H,'flat','HandleVisibility','on');
    H  = findobj(H);
    Hv = get(H,'Visible');
    set(hResAx,'Tag','PermRes','UserData',struct('H',H,'Hv',{Hv}))

    % Draw buttons
    %----------------------------------------------------------------------
    Finter = spm_figure('FindWin','Interactive');
    xyz    = [0;0;0];
    xyz    = spm_XYZreg('RoundCoords',xyz,M,DIM);

    % Create XYZ control objects
    % ---------------------------------------------------------------------
    hFxyz = uicontrol(Finter,'Style','Text',...
            'Position',[010 010 265 030].*WS);
    uicontrol(Finter,'Style','Text','String','co-ordinates',...
            'Position',[020 033 078 016].*WS,...
            'FontAngle','Italic',...
            'FontSize',FS(10),...
            'HorizontalAlignment','Left',...
            'ForegroundColor','w')

    uicontrol(Finter,'Style','Text','String','x =',...
            'Position',[020 015 024 018].*WS,...
            'FontName',PF.times,'FontSize',FS(10),'FontAngle','Italic',...
            'HorizontalAlignment','Center');
    hX   = uicontrol(Finter,'Style','Edit','String',...
            sprintf('%.2f',xyz(1)),...
            'ToolTipString','enter x-coordinate',...
            'Position',[044 015 056 020].*WS,...
            'FontSize',FS(10),'BackGroundColor',[.8,.8,1],...
            'HorizontalAlignment','Right',...
            'Tag','hX',...
            'Callback','spm_bms_display('''',''plot_xyz'')');

    uicontrol(Finter,'Style','Text','String','y =',...
            'Position',[105 015 024 018].*WS,...
            'FontName',PF.times,'FontSize',FS(10),'FontAngle','Italic',...
            'HorizontalAlignment','Center')
    hY   = uicontrol(Finter,'Style','Edit','String',...
            sprintf('%.2f',xyz(2)),...
            'ToolTipString','enter y-coordinate',...
            'Position',[129 015 056 020].*WS,...
            'FontSize',FS(10),'BackGroundColor',[.8,.8,1],...
            'HorizontalAlignment','Right',...
            'Tag','hY',...
            'Callback','spm_bms_display('''',''plot_xyz'')');

    uicontrol(Finter,'Style','Text','String','z =',...
            'Position',[190 015 024 018].*WS,...
            'FontName',PF.times,'FontSize',FS(10),'FontAngle','Italic',...
            'HorizontalAlignment','Center')
    hZ   = uicontrol(Finter,'Style','Edit','String',...
            sprintf('%.2f',xyz(3)),...
            'ToolTipString','enter z-coordinate',...
            'Position',[214 015 056 020].*WS,...
            'FontSize',FS(10),'BackGroundColor',[.8,.8,1],...
            'HorizontalAlignment','Right',...
            'Tag','hZ',...
            'Callback','spm_bms_display('''',''plot_xyz'')');

    % Voxel value reporting pane
    % ---------------------------------------------------------------------
    hFconB = uicontrol(Finter,'Style','Text',...
            'Position',[280 010 110 030].*WS);
    uicontrol(Finter,'Style','Text','String','voxel value',...
            'Position',[285 035 085 016].*WS,...
            'FontAngle','Italic',...
            'FontSize',FS(10),...
            'HorizontalAlignment','Left',...
            'ForegroundColor','w')
    hSPM = uicontrol(Finter,'Style','Text','String','',...
            'Position',[285 012 100 020].*WS,...
            'FontSize',FS(10),...
            'HorizontalAlignment','Center');

    % Store UserData
    % ---------------------------------------------------------------------
    set(hFxyz,'Tag','hFxyz','UserData',struct(...
              'hReg',   [],...
              'M',      M,...
              'iM',     iM,...
              'DIM',    DIM,...
              'XYZ',    xSPM.XYZmm,...
              'Z',      xSPM.Z,...
              'hX',     hX,...
              'hY',     hY,...
              'hZ',     hZ,...
              'hSPM',   hSPM,...
              'xSPM',   xSPM,...
              'fhFxyz', hFxyz,...
              'hMIPax', hMIPax,...
              'xyz',    xyz,...
              'thres',  BMS.xSPM.thres,...
              'scale',  BMS.xSPM.scale,...
              'vols',   BMS.xSPM.vols,...
              'k',      BMS.xSPM.k,...
              'BMS',    BMS));

    set([hX,hY,hZ],'UserData',hFxyz);

    % Register with hReg
    % ---------------------------------------------------------------------
    spm_XYZreg('XReg',hReg,hFxyz,'spm_results_ui');

    % Model partition
    % ---------------------------------------------------------------------
    uicontrol(Finter,'Style','PushButton','String','compare subsets',...
                'FontSize',FS(10),...
                'ToolTipString','Create and compare subsets of models',...
                'Callback','spm_bms_display('''',''partition'');',...
                'Interruptible','on','Enable','on',...
                'Position',[130 055 140 020].*WS,...
                'ForegroundColor','k');

    % Compare groups
    % ---------------------------------------------------------------------
    uicontrol(Finter,'Style','PushButton','String','compare groups',...
                'FontSize',FS(10),...
                'ToolTipString','Compare two groups. E.g. controls vs patients.',...
                'Callback','spm_bms_display('''',''groups'');',...
                'Interruptible','on','Enable','on',...
                'Position',[015 055 100 020].*WS,...
                'ForegroundColor','k');

    % Draw Save, Clear and Exit
    % ---------------------------------------------------------------------
    hClear = uicontrol(Finter,'Style','PushButton','String','clear',...
            'ToolTipString','clears results subpane',...
            'FontSize',FS(9),'ForegroundColor','b',...
            'Callback',['spm_results_ui(''Clear''); ',...
            'spm_input(''!DeleteInputObj''),',...
            'spm_clf(''Satellite'')'],...
            'Interruptible','on','Enable','on',...
            'DeleteFcn','spm_clf(''Graphics'')',...
            'Position',[285 055 035 020].*WS);

    hExit  = uicontrol(Finter,'Style','PushButton','String','exit',...
            'ToolTipString','exit the results section',...
            'FontSize',FS(9),'ForegroundColor','r',...
            'Callback',['spm_clf(''Interactive''),spm_clf(''Graphics''),'...
            'close(spm_figure(''FindWin'',''Satellite'')),'...
            'clear'],...
            'Interruptible','on','Enable','on',...
            'Position',[325 055 035 020].*WS);

    hHelp  = uicontrol(Finter,'Style','PushButton','String','',...
            'ToolTipString','',...
            'FontSize',FS(9),'ForegroundColor','g',...
            'Callback','',...
            'Interruptible','on','Enable','off',...
            'Position',[365 055 020 020].*WS);

    % Change options
    % ---------------------------------------------------------------------
    uicontrol(Finter,'Style','Text',...
            'Position',[125 090 150 085].*WS)
    uicontrol(Finter,'Style','Text','String','options',...
            'Position',[135 168 60 015].*WS,...
            'FontAngle','Italic',...
            'FontSize',FS(10),...
            'HorizontalAlignment','Left',...
            'ForegroundColor','w')
    uicontrol(Finter,'Style','PushButton','String','results',...
            'Position',[130 145 140 020].*WS,...
            'ToolTipString',...
            'BMS Maps (Results)',...
            'Callback','spm_run_bms_vis',...
            'Interruptible','on','Enable','on',...
            'FontSize',FS(10),'ForegroundColor','k')
    uicontrol(Finter,'Style','PushButton','String','change model',...
            'Position',[130 120 140 020].*WS,...
            'ToolTipString',...
            'Change model/data',...
            'Callback','spm_bms_display('''',''change_data'')',...
            'Interruptible','on','Enable','on',...
            'FontSize',FS(10),'ForegroundColor','k')
    uicontrol(Finter,'Style','PushButton','String','threshold',...
            'Position',[130 95 68 020].*WS,...
            'ToolTipString',...
            'Change threshold (same data)',...
            'Callback','spm_bms_display('''',''change_thres'')',...
            'Interruptible','on','Enable','on',...
            'FontSize',FS(8),'ForegroundColor','k')
    uicontrol(Finter,'Style','PushButton','String','scale',...
            'Position',[202 95 68 020].*WS,...
            'ToolTipString',...
            'Change scale (same data)',...
            'Callback','spm_bms_display('''',''change_scale'')',...
            'Interruptible','on','Enable','on',...
            'FontSize',FS(8),'ForegroundColor','k')

    %-p-values
    %------------------------------------------------------------------
    uicontrol(Finter,'Style','Text','String','p-values',...
            'Position',[020 168 050 015].*WS,...
            'FontAngle','Italic',...
            'FontSize',FS(10),...
            'HorizontalAlignment','Left',...
            'ForegroundColor','w')
    uicontrol(Finter,'Style','PushButton','String','whole brain','FontSize',FS(10),...
            'ToolTipString',...
            'tabulate summary of local maxima, p-values & statistics',...
            'Callback','spm_bms_display('''',''list'');',...
            'Interruptible','on','Enable','on',...
            'Position',[015 145 100 020].*WS)
    uicontrol(Finter,'Style','PushButton','String','current cluster','FontSize',FS(10),...
            'ToolTipString',...
            'tabulate p-values & statistics for local maxima of nearest cluster',...
            'Callback','spm_bms_display('''',''listCluster'');',...
            'Interruptible','on','Enable','on',...
            'Position',[015 120 100 020].*WS)
    uicontrol(Finter,'Style','PushButton','String','ROI','FontSize',FS(10),...
            'ToolTipString',...
            'plot probabilities for selected ROI',...
            'Callback','spm_bms_display('''',''plotROI'');',...
            'Position',[015 095 100 020].*WS)

    uicontrol(Finter,'Style','PushButton','String','small volume','FontSize',FS(10),...
            'ToolTipString',['Small Volume Correction - probability values ',...
            'for a small search region'],...
            'Callback','spm_bms_display('''',''listVOI'');',...
            'Interruptible','on','Enable','on',...
            'Position',[015 095 100 020].*WS)

    % Draw Options
    % ---------------------------------------------------------------------
uicontrol(Finter,'Style','Text',...
            'Position',[280 090 110 085].*WS)
    uicontrol(Finter,'Style','Text','String','display',...
            'Position',[290 168 065 015].*WS,...
            'FontAngle','Italic',...
            'FontSize',FS(10),...
            'HorizontalAlignment','Left',...
            'ForegroundColor','w')

    strp  = { 'plot...','current voxel','ROI'};
    tstrp = { 'plot results from comparisons at: ',...
        'current voxel / ','region of interest /'};
    tmpp  = { 'spm_bms_display('''',''do_plot'')',...
            'spm_bms_display('''',''do_ROI'')'};

    uicontrol(Finter,'Style','PopUp','String',strp,'FontSize',FS(10),...
            'ToolTipString',cat(2,tstrp{:}),...
            'UserData',tmpp,...
            'Callback','spm_bms_display('''',''overlays'')',...
            'Interruptible','on','Enable','on',...
            'Position',[285 145 100 020].*WS)
    str  = { 'overlays...','slices','sections','render','previous sections'};
    tstr = { 'overlay results on another image: ',...
        '3 slices / ''ortho sections / ','render /','previous ortho sections'};
    tmp  = { 'spm_transverse(''set'',xSPM,hReg)',...
            'spm_sections(xSPM,hReg);global st;st.vols{1}.blobs{1}.min=xSPM.u;spm_orthviews(''redraw'');',...
            ['spm_render(   struct( ''XYZ'',    xSPM.XYZ,',...
            '''t'',     xSPM.Z'',',...
            '''mat'',   xSPM.M,',...
            '''dim'',   xSPM.DIM))'],...
            ['global prevsect;','spm_sections(xSPM,hReg,prevsect)'],...
            ['global prevrend;','if ~isstruct(prevrend)',...
            'prevrend = struct(''rendfile'','''',''brt'',[],''col'',[]); end;',...
            'spm_render(    struct( ''XYZ'',    xSPM.XYZ,',...
            '''t'',     xSPM.Z'',',...
            '''mat'',   xSPM.M,',...
            '''dim'',   xSPM.DIM),prevrend.brt,prevrend.rendfile)']};
    uicontrol(Finter,'Style','PopUp','String',str,'FontSize',FS(10),...
            'ToolTipString',cat(2,tstr{:}),...
            'Callback','spm_bms_display('''',''overlays'')',...
            'UserData',tmp,...
            'Interruptible','on','Enable','on',...
            'Position',[285 120 100 020].*WS)
    uicontrol(Finter,'Style','PushButton','String','save','FontSize',...
            FS(10),'ToolTipString','save thresholded BMS as image',...
            'Callback','spm_bms_display('''',''save'')',...
            'Interruptible','on','Enable','on',...
            'Position',[285 095 100 020].*WS)
    user_data = get(hFxyz,'UserData');
    set(Finter,'UserData',user_data,...
            'HandleVisibility','callback')

    varargout = { hReg };

    % Do plot - action: 'do_plot'
    % =====================================================================
    case 'do_plot'

        fig       = gcf;
        user_data = get(fig,'UserData');
        iM        = user_data.iM;
        BMS       = user_data.BMS;
        user_data = get(user_data.hMIPax,'UserData');
        user_data = get(user_data.hMIPxyz,'UserData');
        xyz_vx    = iM*[user_data; 1];
        spm_bms_display_vox(BMS,xyz_vx(1:3));

    % Do ROI - action: 'do_ROI'
    % =====================================================================
    case 'do_ROI'

        fig       = gcf;
        user_data = get(fig,'UserData');
        iM        = user_data.iM;
        BMS       = user_data.BMS;
        user_data = get(user_data.hMIPax,'UserData');
        user_data = get(user_data.hMIPxyz,'UserData');
        xyz_vx    = iM*[user_data; 1];
        spm_bms_display_ROI(BMS);

    % Do overlays - action: 'overlays'
    % =====================================================================
    case 'overlays'

        h   = gcbo;
        v   = get(h,'Value');
        if v==1, return, end
        set(h,'Value',1)
        CBs       = get(h,'UserData');
        fig       = gcf;
        user_data = get(fig,'UserData');
        hReg      = user_data.hReg;
        xSPM      = user_data.xSPM;
        eval(CBs{v-1})

    % Save - action: 'save'
    % =====================================================================
    case 'save'

        fig       = gcf;
        user_data = get(fig,'UserData');
        xSPM      = user_data.xSPM;
        spm_write_filtered(xSPM.Z,xSPM.XYZ,xSPM.DIM,xSPM.M,'Results saved');

    % Plot xyz BMS values - action: 'plot_xyz'
    % =====================================================================
    case 'plot_xyz'

        hC    = gcbo;
        d     = find(strcmp(get(hC,'Tag'),{'hX','hY','hZ'}));
        hFxyz = get(hC,'UserData');
        UD    = get(hFxyz,'UserData');
        xyz   = UD.xyz;
        nxyz  = xyz;

        o = evalin('base',['[',get(hC,'String'),']'],'sprintf(''error'')');
        if ischar(o) || length(o)>1
            warning('%s: Error evaluating ordinate:\n\t%s',...
                mfilename,lasterr)
        else
            nxyz(d) = o;
            nxyz = spm_XYZreg('RoundCoords',nxyz,UD.M,UD.DIM);
        end

        if abs(xyz(d)-nxyz(d))>0
            UD.xyz = nxyz; set(hFxyz,'UserData',UD)
            if ~isempty(UD.hReg), spm_XYZreg('SetCoords',nxyz,UD.hReg,hFxyz); end
            set(hC,'String',sprintf('%.3f',nxyz(d)))
            i  = spm_XYZreg('FindXYZ',UD.xyz,UD.XYZ);
            if isempty(i), str = ''; else str = sprintf('%6.2f',UD.Z(i)); end
               set(UD.hSPM,'String',str);
        end

        fig = gcf;
        set(fig,'UserData',UD)

    % Change model
    % =====================================================================
    case 'change_data'

        fig         = gcf;
        user_data   = get(fig,'UserData');
        job.img{1}  = '';
        job.file{1} = user_data.BMS.fname;
        job.thres   = user_data.thres;
        job.scale   = user_data.scale;
        job.k       = user_data.k;
        spm_run_bms_vis(job);

    % Change threshold
    % =====================================================================
    case 'change_thres'

        fig         = gcf;
        user_data   = get(fig,'UserData');
        job.img{1}  = user_data.vols;
        job.thres   = [];
        job.file{1} = user_data.BMS.fname;
        job.scale   = user_data.scale;
        job.k       = user_data.k;
        spm_run_bms_vis(job);

    % Change scale
    % =====================================================================
    case 'change_scale'

        fig         = gcf;
        user_data   = get(fig,'UserData');
        job.img{1}  = user_data.vols;
        job.thres   = user_data.thres;
        job.scale   = [];
        job.k       = user_data.k;
        job.file{1} = user_data.BMS.fname;
        spm_run_bms_vis(job);

    % List p-values
    % =====================================================================
    case 'list'

        fig       = gcf;
        user_data = get(fig,'UserData');
        xSPM      = user_data.xSPM;
        hReg      = user_data.hReg;
        xSPM.STAT = 'P';
        xSPM.Z    = xSPM.z_ps;
        spm_list('List',xSPM,hReg);

    % List cluster
    % =====================================================================
    case 'listCluster'

        fig       = gcf;
        user_data = get(fig,'UserData');
        xSPM      = user_data.xSPM;
        hReg      = user_data.hReg;
        xSPM.STAT = 'P';
        xSPM.Z    = xSPM.z_ps;
        spm_list('ListCluster',xSPM,hReg);

    % Small volume
    % =====================================================================
    case 'listVOI'

        fig       = gcf;
        user_data = get(fig,'UserData');
        xSPM      = user_data.xSPM;
        hReg      = user_data.hReg;
        xSPM.STAT = 'P';
        spm_VOI(xSPM.SPM,xSPM,hReg);

    % Model partitioning
    % =====================================================================
    case 'partition'
        fig       = gcf;
        user_data = get(fig,'UserData');
        BMS       = user_data.BMS;
        spm_bms_partition(BMS);

    % Model partitioning
    % =====================================================================
    case 'groups'

        fig         = gcf;
        user_data   = get(fig,'UserData');
        BMS         = user_data.BMS;

        con_image   = spm_bms_compare_groups;

        job.img{1}  = con_image;
        job.file{1} = user_data.BMS.fname;
        job.thres   = [];
        job.scale   = [];
        job.k       = [];
        spm_run_bms_vis(job);

end  % End switch

end  % End function

author: spm


text:
function [edges,idx,edgemap]=uniqedges(elem)
%
% [edges,idx,edgemap]=uniqedges(elem)
%
% return the unique edge list from a surface or tetrahedral mesh
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
%
% input:
%     elem: a list of elements, each row is a list of nodes for an element.
%           elem can have 2, 3 or 4 columns
%
% output:
%     edge: unique edges in the mesh, denoted by a pair of node indices
%     idx:  index of the output in the raw edge list (returned by meshedge)
%     edgemap: index of the raw edges in the output list (for triangular mesh)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(size(elem)==2)
   edges=elem;
elseif(size(elem)>=3)
   edges=meshedge(elem);
else
   error('invalid input');
end

[uedges,idx,jdx]=unique(sort(edges,2),'rows');
edges=edges(idx,:);
edgemap=reshape(jdx,size(elem));

repo: REST


text:
function z = geoatan2(x, y)
%-------------------------------------------------------
% z = geoatan2(x, y);
% Yudan Yi, May 26, 2005
%-------------------------------------------------------
% return -pi~pi of y/x
z = 0.0;
if (abs(x)<1e-12)  if (y>0) z = pi/2.0; else z =-pi/2.0; end; return; end;
z = atan(abs(y/x));
if (x>0)
	if (y<0) z =-z; end; return;
else
	if (y>0) z = pi-z; else z =-pi+z; end;
end
return;

repo: ekfmonoslam


text:
function [Btuph] = GW2Btuph(GW)
% Convert power from gigawatts to British thermal units per hour.
% Chad A. Greene 2012
Btuph = GW*3412141633.1279;


text:
function contrib_prev = reverse_max_vl(l_prev, contrib_curr)
    contrib_prev = vl_nnpool(l_prev, [2,2], contrib_curr, 'Stride', 2);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function Rank=ks(X)
%+++ Employ the K-S algorithm for selecting the representative samples;
%+++ Basic idea of KS: Find from the candidates the sample whose minimal
%    distance from the samples in representative sets is the maximal.
%+++ Hongdong Li, May 10,2008.
%+++ Revised in Dec. 5, 2009.

tic;
[Mx,Nx]=size(X);
Rank=zeros(1,Mx);           %+++ Record the indices of the representative samples.
out=1:Mx;
D=distli(X);
[i j]=find(D==max(max(D)));
Rank(1)=i(1);               %+++ Initializes as the two samples of the furthest distance.
Rank(2)=j(1);
out([i(1) j(1)])=[];        %+++ The remaining samples.
%+++ Iteration of  K-S algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iter=3;
while iter<=Mx
   in=Rank(find(Rank>0));
   Dsub=D(in,out);
   [minD,indexmin]=min(Dsub);
   [maxD,indexmax]=max(minD);
   Vadd=out(indexmax);
   Rank(iter)=Vadd;
   out(find(out==Vadd))=[];
   iter=iter+1;
end
toc;
%+++ Iteration ended %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%% END OF SUB
function D=distli(X)
X=X';
[D,N] = size(X);
X2 = sum(X.^2,1);
D = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;

author: viggin

repo: domain-adaptation-toolbox