text:
function reslicedImage = reslice(this, targetGeometry, varargin)
% Reslices image to image geometry of other image using spm_reslice
%   Y = MrImageSpm4D()
%   resliceImage = Y.reslice(targetGeometry, ...
%       'spmParameterName1', spmParameterValue1, ...
%       ...
%       'spmParameterNameN', spmParameterValueN)
%
%   OR
%   Y.reslice(otherImage);
%
% This is a method of class MrImageSpm4D.
%
% IN
%   targetGeometry  object of MrImageGeometry or MrImage
%                   Image will be resliced to this geometry
%   interpolation   degree of b-spline interpolation for estimation and reslicing
%                   default: 7
%   wrapping        fold-over direction (phase encode)
%                   default: [0 0 0] % none
%   masking         mask incomplete timeseries?
%                   default: true
%
%
% OUT
%
%   reslicedImage
%
% EXAMPLE
%   Y = MrImage();
%   Z = MrImage();
%   targetGeometry = Z.geometry;
%   Y.reslice(targetGeometry)
%
%   See also MrImage MrImageGeometry spm_reslice spm_run_coreg

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-18
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

defaults.interpolation = 7;     % degree of b-spline interpolation for estimation and reslicing
defaults.wrapping = [0 0 0];    % fold-over direction (phase encode)
defaults.masking = 1;           % mask incomplete timeseries?

args = tapas_uniqc_propval(varargin, defaults);
reslicedImage = this.copyobj();

% Save as nifti to use spm functionality
% but check if file already exists, give new filename to prevent
% accidental overwrite
if isnumeric(reslicedImage.parameters.save.keepCreatedFiles)
    keepCreatedFiles = reslicedImage.parameters.save.keepCreatedFiles;
else
    keepCreatedFiles = ~strcmpi(reslicedImage.parameters.save.keepCreatedFiles, 'none');
end

changeFilename = isfile(reslicedImage.get_filename) && ~keepCreatedFiles;
if changeFilename
    origFilename = reslicedImage.parameters.save.fileName;
    [~, tmpName] = fileparts(tempname);
    reslicedImage.parameters.save.fileName = [tmpName, '.nii'];
end

reslicedImage.save('fileName', reslicedImage.get_filename('prefix', 'raw'));

% check whether input is actually a geometry
isGeometry = isa(targetGeometry, 'MrImageGeometry');
if ~isGeometry
    if isa(targetGeometry, 'MrImage')
        targetGeometry = targetGeometry.geometry;
    else
        disp('Input has to be of class MrImage or MrImageGeometry.');
    end
end

[~, ~, isEqualGeom3D] = targetGeometry.diffobj(reslicedImage.geometry);

if ~isEqualGeom3D

    % Dummy 3D image with right geometry is needed for resizing
    emptyImage = targetGeometry.create_empty_image('z', 1);
    emptyImage.parameters.save.path = reslicedImage.parameters.save.path;
    emptyImage.save();
    fnTargetGeometryImage = emptyImage.get_filename;

    matlabbatch = reslicedImage.get_matlabbatch('reslice', fnTargetGeometryImage, args);
    save(fullfile(reslicedImage.parameters.save.path, 'matlabbatch.mat'), ...
        'matlabbatch');
    spm_jobman('run', matlabbatch);

    % clean up: move/delete processed spm files, load new data into matrix
    reslicedImage.finish_processing_step('reslice', fnTargetGeometryImage);
end
% set back to original filename
if changeFilename
    reslicedImage.parameters.save.fileName = origFilename;
end

end

author: translationalneuromodeling

repo: tapas


function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu,ww)
%function var=rfvar3(ydata,lags,xdata,breaks,lambda,mu)
% This algorithm goes for accuracy without worrying about memory requirements.
% ydata:   dependent variable data matrix
% xdata:   exogenous variable data matrix
% lags:    number of lags
% breaks:  rows in ydata and xdata after which there is a break.  This allows for
%          discontinuities in the data (e.g. war years) and for the possibility of
%          adding dummy observations to implement a prior.  This must be a column vector.
%          Note that a single dummy observation becomes lags+1 rows of the data matrix,
%          with a break separating it from the rest of the data.  The function treats the
%          first lags observations at the top and after each "break" in ydata and xdata as
%          initial conditions.
% lambda:  weight on "co-persistence" prior dummy observations.  This expresses
%          belief that when data on *all* y's are stable at their initial levels, they will
%          tend to persist at that level.  lambda=5 is a reasonable first try.  With lambda<0,
%          constant term is not included in the dummy observation, so that stationary models
%          with means equal to initial ybar do not fit the prior mean.  With lambda>0, the prior
%          implies that large constants are unlikely if unit roots are present.
% mu:      weight on "own persistence" prior dummy observation.  Expresses belief
%          that when y_i has been stable at its initial level, it will tend to persist
%          at that level, regardless of the values of other variables.  There is
%          one of these for each variable.  A reasonable first guess is mu=2.
%      The program assumes that the first lags rows of ydata and xdata are real data, not dummies.
%      Dummy observations should go at the end, if any.  If pre-sample x's are not available,
%      repeating the initial xdata(lags+1,:) row or copying xdata(lags+1:2*lags,:) into
%      xdata(1:lags,:) are reasonable subsititutes.  These values are used in forming the
%      persistence priors.

% Original file downloaded from:
% http://sims.princeton.edu/yftp/VARtools/matlab/rfvar3.m

if nargin<7
   scale_ = 0;
else
    % correct for heteroskedasticity
    scale_=1;
end

[T,nvar] = size(ydata);
nox = isempty(xdata);
if ~nox
    [T2,nx] = size(xdata);
else
    T2 = T;
    nx = 0;
    xdata = zeros(T2,0);
end
% note that x must be same length as y, even though first part of x will not be used.
% This is so that the lags parameter can be changed without reshaping the xdata matrix.
if T2 ~= T, error('Mismatch of x and y data lengths'),end
if nargin < 4
    nbreaks = 0;
    breaks = [];
else
    nbreaks = length(breaks);
end
breaks = [0;breaks;T];
smpl = [];
for nb = 1:nbreaks+1
    smpl = [smpl;[breaks(nb)+lags+1:breaks(nb+1)]'];
end
Tsmpl = size(smpl,1);
X = zeros(Tsmpl,nvar,lags);
for is = 1:length(smpl)
    X(is,:,:) = ydata(smpl(is)-(1:lags),:)';
end
X = [X(:,:) xdata(smpl,:)];
y = ydata(smpl,:);

% rescale if heteroskedasticity corrected
if scale_ ==1
    ww = [ww; ones(length(y)-length(ww),1) ];
    y = y./ repmat(ww,1,size(y,2)) ;
    X = X ./ repmat(ww,1,size(X,2));
end
% Everything now set up with input data for y=Xb+e

% Add persistence dummies
if lambda ~= 0 || mu > 0
    ybar = mean(ydata(1:lags,:),1);
    if ~nox
        xbar = mean(xdata(1:lags,:),1);
    else
        xbar = [];
    end
    if lambda ~= 0
        if lambda>0
            xdum = lambda*[repmat(ybar,1,lags) xbar];
        else
            lambda = -lambda;
            xdum = lambda*[repmat(ybar,1,lags) zeros(size(xbar))];
        end
        ydum = zeros(1,nvar);
        ydum(1,:) = lambda*ybar;
        y = [y;ydum];
        X = [X;xdum];
    end
    if mu>0
        xdum = [repmat(diag(ybar),1,lags) zeros(nvar,nx)]*mu;
        ydum = mu*diag(ybar);
        X = [X;xdum];
        y = [y;ydum];
    end
end

% Compute OLS regression and residuals
[vl,d,vr] = svd(X,0);
di = 1./diag(d);
B = (vr.*repmat(di',nvar*lags+nx,1))*vl'*y;
u = y-X*B;
xxi = vr.*repmat(di',nvar*lags+nx,1);
xxi = xxi*xxi';

var.B = B;
var.u = u;
var.xxi = xxi;
var.y   = y;
var.X   = X;


text:
function [u,p,info] = mgDarcy(M,B,f,g,elem,option,varargin)
%% MGDARCY multigrid solvers for Darcy system discretized by RT0-P0 element
%
% [u,p,info] = mgDarcy(M,B,f,g,elem) solves saddle point problem
% discretized from RT0-P0 mixed FEM for Darcy equation.
%
%      |M  B'| |u|  = |f|
%      |B  0 | |p|  = |g|
%
%  A V-cycle multigrid using overlapping Schwarz smoother is implemented.
%  In the first step, mgdivDarcy is called to find initial u satisfying Bu
%  = g. Then at each vertex patch, a local problem with prescribed boundary
%  flux is solved. Detailed description and convergence analysis can be
%  found in the reference.
%
%  It is around 10 times slower than tripremixPoisson since a large for
%  loop is not efficient in MATLAB. In operation count, this is superior.
%  See the complexity analysis in page 18 of the reference.
%
% Reference: L. Chen. Multigrid Methods for Constrained Minimization
% Problems and Application to Saddle Point Problems
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

time = cputime;
%% Size of systems
Nu = length(f);
Np = length(g);
N = max(elem(:));                  % number of nodes
% NT = size(elem,1);               % number of elements

%% Options
% Assign default values to unspecified parameters
if ~exist('option','var'),
    option = [];
end
option = mgoptions(option,Nu);    % parameters
u0 = option.x0;
N0 = option.N0;
tol = option.tol;
maxIt = option.solvermaxit;
mu = option.smoothingstep;
coarsegridsolver = option.coarsegridsolver;
printlevel = option.printlevel;
freeEdge = option.freeEdge;

%% Hierarchical Structure of Mesh
HB = zeros(N,3);
level = 20;
NL(level+1) = N; % now NL(1:level) = 0;
elemi = cell(level,1);
elemi{level} = elem;
freeEdgei = cell(level,1);
freeEdgei{level} = freeEdge;
Pro_u = cell(level,1);
Pro_p = cell(level,1);
for k = level: -1 : 2
    switch option.refType
        case 'red'
            [elemi{k-1},newHB] = uniformcoarsenred(elemi{k}); % coasen red refinement
            if ~isempty(newHB)
                [Pro_u{k-1}, freeEdgei{k-1}] = transferedgered(elemi{k-1},elemi{k},freeEdgei{k}); % transfer operator of u
                Pro_p{k-1} = repmat(speye(size(elemi{k-1},1)),4,1);   % transfer operator of p
            end
        case 'bisect'
            % merge two coarsen of bisection grids s.t. the ratio is 1/4.
            [tempelem,newHB1,tree] = uniformcoarsen(elemi{k}); % coarse bisection
            if ~isempty(newHB1)
                % first coarsen
                [tempPro_u,tempFreeEdge] = transferedgecoarsen(tempelem,elemi{k},tree,freeEdgei{k});
                tempPro_p = transferelem(tempelem,elemi{k},tree);
                % second coarsen
                [elemi{k-1},newHB2,tree] = uniformcoarsen(tempelem); % coarse bisection
                [Pro_u{k-1}, freeEdgei{k-1}] = transferedgecoarsen(elemi{k-1},tempelem,tree,tempFreeEdge);
                Pro_u{k-1} = tempPro_u*Pro_u{k-1};
                Pro_p{k-1} = transferelem(elemi{k-1},tempelem,tree);
                Pro_p{k-1} = tempPro_p*Pro_p{k-1};
                newHB = [newHB1; newHB2];
            end
    end
    if (isempty(newHB)) || (size(elemi{k},1)< 2*N0)
        % no nodes are removed or it reaches the coarsest level
        NL = NL(k:end);
        break;
    end
    NL(k) = NL(k+1) - size(newHB,1); % update NL(k)
    HB(NL(k)+1:NL(k+1),1:3) = newHB(:,1:3);
end
level = length(NL)-1;    % actual level
elemi = elemi(end-level+1: end);
Pro_u = Pro_u(end-level+1: end);
Pro_p = Pro_p(end-level+1: end);
freeEdgei = freeEdgei(end-level+1: end);
% generate edge, elem2edge and freeNode etc
edgei = cell(level,1);
elem2edgei = cell(level,1);
for k = 1:level
    [elem2edgei{k},edgei{k}] = dofedge(elemi{k});
end

%% Matrices in each level
oldf = f;
Mi = cell(level,1);
Bi = cell(level,1);
if size(M,1) > length(freeEdge) % truncate to free edge only
    Mi{level} = M(freeEdge,freeEdge);
    Bi{level} = B(:,freeEdge);
    f  = f(freeEdge);
    u0 = u0(freeEdge);
else
    Mi{level} = M;
    Bi{level} = B;
end
Res_u = cell(level,1);
Res_p = cell(level,1);
for k = level:-1:2
    Res_u{k} = transpose(Pro_u{k-1});
    Res_p{k} = transpose(Pro_p{k-1});
    Mi{k-1} = Res_u{k}*Mi{k}*Pro_u{k-1};
    Bi{k-1} = Res_p{k}*Bi{k}*Pro_u{k-1};
end

% %% Exact solver: for debug only
% C = sparse(size(B,1),size(B,1));
% A = [M B';B C];
% F = [f; zeros(size(B,1),1)];
% tempu = A\F;
% exactSigma = tempu(1:Ndof);

%% MG cycles
% initial set up
k = 1;
f0 = f - Mi{level}*u0;
g0 = g - Bi{level}*u0;
% find u satisfy Bu = g
u = u0 + mgdivDarcy(f0,g0);
% temp = abs(g-Bi{level}*u);
% idx = temp<1e-14;
% findelem(node,elem,idx,'noindex','FaceColor','c');
if printlevel >= 1
    fprintf('Multigrid Vcycle Iteration \n')
end
err = zeros(maxIt,1);
err(1) = 1;
r0 = f - Mi{level}*u;
r = r0;
while (max(err(k,:)) > tol) && (k <= maxIt)
    k = k + 1;
    % Step 2: Compute Br by one Vcylce MG
    Br = vcycle(r);
    % Step 3: Correct the solution
    u = u + Br;
    % Step 1: Form residual r
    r = r - Mi{level}*Br;
    err(k) = sqrt(abs(Br'*r/(u'*r0))); % approximate relative error in energy norm
    if printlevel >= 2
        fprintf('#dof: %8.0u,  #nnz: %8.0u, MG Vcycle iter: %2.0u, err = %8.4e\n',...
                 Nu+N, nnz(M)+nnz(B), k-1, err(k));
    end
end

err = err(1:k,:);
itStep = k-1;

%% Find pressure
Ap = B*B';
ubar = zeros(Nu,1);
ubar(freeEdge) = u;
b = B*(oldf-M*ubar);
p = zeros(Np,1);
if option.isPureNeumannBC
    freep = 1:Np-1;
else
    freep = 1:Np;
end
p(freep) = Ap(freep,freep)\b(freep);

%% Output
time = cputime - time;
if k > maxIt
    flag = 1;
else
    flag = 0;
end
if printlevel >= 2
    fprintf('#dof: %8.0u, level: %2.0u,   coarse grid %2.0u, #nnz: %8.0u\n',...
              Nu+N, level, size(Mi{1},1), nnz(Mi{1}))
end
if printlevel >= 1
    fprintf('#dof: %8.0u,  #nnz: %8.0u, smoothing: %2.0u, iter: %2.0u,   err = %8.4e,   time = %4.2g s\n',...
                 Nu+Np, nnz(M)+2*nnz(B), mu, itStep, err(end), time)
end
if (flag == 1) && (printlevel>0)
   fprintf('NOTE: the iterative method does not converge! \n');
end
info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions vcycle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Vcycle MG
    function Br = vcycle(r,J)        % solve equations Ae = r in each level
    if nargin<=1
        J = level;
    end
    ri = cell(J,1);            % record residual in each level
    ei = cell(J,1);            % record err in each level
    ri{J} = r;
    for i = J:-1:2
        ei{i} = zeros(size(Mi{i},1),1);
        ei{i} = SchwarzsmootherDarcy(Mi{i},Bi{i},ri{i},zeros(size(Bi{i},1),1),...
                ei{i},elemi{i},edgei{i},elem2edgei{i},freeEdgei{i},mu);
        ri{i-1} = Res_u{i}*(ri{i} - Mi{i}*ei{i});
    end
    if strcmp(coarsegridsolver,'direct')
        Ncoarse = size(Bi{1},1);
        C1 = sparse(Ncoarse,Ncoarse);
        A1 = [Mi{1} Bi{1}'; Bi{1} C1];
        F1 = [ri{1}; zeros(Ncoarse,1)];
        if option.isPureNeumannBC
            bigu = A1(1:end-1,1:end-1)\F1(1:end-1); % pure Neumann
        else
            bigu = A1\F1;
        end
        ei{1} = bigu(1:size(Mi{1},1));
    end
    for i = 2:J
        ei{i} = ei{i} + Pro_u{i-1}*ei{i-1};
        ei{i} = SchwarzsmootherDarcy(Mi{i},Bi{i},ri{i},zeros(size(Bi{i},1),1),...
                ei{i},elemi{i},edgei{i},elem2edgei{i},freeEdgei{i},-mu);
    end
    Br = ei{J};
    end
%% div Darcy
% Use one V-cycle with post-smoothing only to find u s.t. div u = g

    function u = mgdivDarcy(rf,rg)
    J = level;
    rfi = cell(J,1);
    rgi = cell(J,1);
    eui = cell(J,1);
    rfi{J} = rf;
    rgi{J} = rg;
    % restrict the residual to the coarse level
    for i = J:-1:2
        rfi{i-1} = Res_u{i}*rfi{i};
        rgi{i-1} = Res_p{i}*rgi{i};
    end
    % exact solve in the coarest mesh
    Ncoarse = size(Bi{1},1);
    C1 = sparse(Ncoarse,Ncoarse);
    A1 = [Mi{1} Bi{1}'; Bi{1} C1];
    F1 = [rfi{1}; rgi{1}];
    if option.isPureNeumannBC
        bigu = A1(1:end-1,1:end-1)\F1(1:end-1); % pure Neumann
    else
        bigu = A1\F1;
    end
    eui{1} = bigu(1:size(Mi{1},1));
    % prolongate to the fine level and post-smoothing in each element
    for i = 2:J
        eui{i} = Pro_u{i-1}*eui{i-1};
        eui{i} = SchwarzsmootherelemDarcy(Mi{i},Bi{i},rfi{i},rgi{i},...
                eui{i},elemi{i},edgei{i},elem2edgei{i},freeEdgei{i},1);
    end
    u = eui{J};
    end
end

repo: ifem


text:
run vlfeat/toolbox/vl_setup ;
run matconvnet/matlab/vl_setupnn ;
addpath('core');
addpath('helpers');

author: aravindhm

repo: deep-goggle


text:
function vn = backward_chain(w,tr,tc,la)
%
%Backward chain from goal of given world to find the initial states that
%are reachable
%
t.r = tr;
t.c = tc;
t.l = la;
q = qu_new;% start node queue
q = qu_enqu(q,t);% enqueue goal
v(1:size(w,1),1:size(w,2),1:la) = 0;
while qu_empty(q) ~= 1
    cu = qu_front(q);
    q = qu_dequ(q);% dequeue goal
    if (w(cu.r,cu.c,cu.l) ~= 1)&&(cu.l~=1)% check if obstacle or initial state
        if v(cu.r,cu.c,cu.l-1) ~= 1
            t.r = cu.r;
            t.c = cu.c;
            t.l = cu.l - 1;
            q = qu_enqu(q,t); % enqueue descendant
            v(cu.r,cu.c,cu.l-1) = 1;% mark descendant as visited
        end
        if (cu.r > 1)&&(v(cu.r-1,cu.c,cu.l-1) ~= 1)
            t.r = cu.r - 1;
            t.c = cu.c;
            t.l = cu.l - 1;
            q = qu_enqu(q,t);
            v(cu.r-1,cu.c,cu.l-1) = 1;
        end
        if (cu.c > 1)&&(v(cu.r,cu.c-1,cu.l-1) ~= 1)
            t.r = cu.r;
            t.c = cu.c - 1;
            t.l = cu.l - 1;
            q = qu_enqu(q,t);
            v(cu.r,cu.c-1,cu.l-1) = 1;
        end
        if (cu.r < size(w,1))&&(v(cu.r+1,cu.c,cu.l-1) ~= 1)
            t.r = cu.r + 1;
            t.c = cu.c;
            t.l = cu.l - 1;
            q = qu_enqu(q,t);
            v(cu.r+1,cu.c,cu.l-1) = 1;
        end
        if (cu.c < size(w,2))&&(v(cu.r,cu.c+1,cu.l-1) ~= 1)
            t.r = cu.r;
            t.c = cu.c + 1;
            t.l = cu.l - 1;
            q = qu_enqu(q,t);
            v(cu.r,cu.c+1,cu.l-1) = 1;
        end
    end
end
q = qu_free(q);% free queue space
vn = v; % return values
end

author: Sable

repo: mcbench-benchmarks


text:
function [ past_tracks, future_tracks ] = get_tracks( curr_tracks, curr_tags, curr_date_index )
%GET_TRACKS Summary of this function goes here
%   Detailed explanation goes here

% Getting the current eddy indexes in tracks and track length
track_lengths = cellfun(@(x) size(x, 1), curr_tracks);

curr_eddy_indexes = zeros(size(track_lengths));
for i = 1:length(curr_tracks)
    curr_eddy_indexes(i) = find(curr_tracks{i}(:, 3) == curr_date_index, 1);
end

% Pre-allocate past_tracks and future_tracks
past_tracks = nan(sum(curr_eddy_indexes) + length(curr_tracks), 3);
future_tracks = nan(sum(track_lengths) - sum(curr_eddy_indexes) +  2 * length(curr_tracks), 3);

next_future_track_index = 1;
next_past_track_index = 1;
for i = 1:length(curr_tracks)
    curr_track = curr_tracks{i};
    curr_tag = curr_tags{i};
    curr_eddy_index = curr_eddy_indexes(i);

    past_tracks(next_past_track_index:(next_past_track_index + curr_eddy_index - 1), 1:2) = ...
        curr_track(1:curr_eddy_index, 1:2);
    past_tracks(next_past_track_index:(next_past_track_index + curr_eddy_index - 1), 3) = ...
        curr_tag(1:curr_eddy_index);
    next_past_track_index = next_past_track_index + curr_eddy_index + 1;

    future_tracks(next_future_track_index:(next_future_track_index + track_lengths(i) - curr_eddy_index), 1:2) = ...
        curr_track(curr_eddy_index:end, 1:2);

    future_tracks(next_future_track_index:(next_future_track_index + track_lengths(i) - curr_eddy_index), 3) = ...
        curr_tag(curr_eddy_index:end);
    next_future_track_index = next_future_track_index + track_lengths(i) - curr_eddy_index + 2;
end

end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
% matperm() - transpose and sign rows of x to match y (run after matcorr() )
%
% Usage: >> [permx indperm] = matperm(x,y,indx,indy,corr);
%
% Inputs:
%   x     = first input matrix
%   y     = matrix with same number of columns as x
%   indx  = column containing row indices for x (from matcorr())
%   indy  = column containing row indices for y (from matcorr())
%   corr  = column of correlations between indexed rows of x,y (from matcorr())
%           (used only for its signs, +/-)
% Outputs:
%   permx   = the matrix x permuted and signed according to (indx, indy,corr)
%             to best match y. Rows of 0s added to x to match size of y if nec.
%   indperm = permutation index turning x into y;
%
% Authors: Scott Makeig, Sigurd Enghoff & Tzyy-Ping Jung
%          SCCN/INC/UCSD, La Jolla, 2000

% Copyright (C) 1996 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 04-22-99  Adjusted for fixes and speed by Sigurd Enghoff & Tzyy-Ping Jung
% 01-25-02  Reformated help & license, added links -ad

function [permx,indperm]= matperm(x,y,indx,indy,corr)

[m,n] = size(x);
[p,q] = size(y);
[ix,z] = size(indx);
[iy,z] = size(indy);
oldm = m;

errcode=0;
if  ix ~= iy | p ~= iy,
 fprintf('matperm: indx and indy must be column vectors, same height as y.\n');
 errcode=1
end;

if n~=q,
   fprintf('matperm(): two matrices must be same number of columns.\n');
   errcode=2;
else
  if m<p,
  		x = [x;zeros(p-m,n)];	% add rows to x to match height of y
  		p=m;
  elseif p<m,
  		y = [y;zeros(m-p,n)];	% add rows to y to match height of x
  		m=p;
  end;
end;
if errcode==0,
%
% Return the row permutation of matrix x most correlated with matrix y:
%  plus the resulting permutation index
%
  indperm = [1:length(indx)]';	% column vector [1 2 ...nrows]
  permx  = x(indx,:);
  indperm = indperm(indx,:);
  ydni(indy) = 1:length(indy);
  permx = permx(ydni,:);% put x in y row-order
  indperm = indperm(ydni,:);
  permx = permx.*(sgn(corr(ydni))*ones(1,size(permx,2)));
								 % make x signs agree with y
  permx = permx(1:oldm,:);		 % throw out bottom rows if
								 % they were added to match y
  indperm = indperm(1:oldm,:);
end;

return

function vals=sgn(data)

 vals = 2*(data>=0)-1;

return

author: PatternRecognition

repo: OpenBMI


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
function [sys,x0,str,ts] = quadrotor_plot(t,x,u,flag,s,plot,enable,vehicle)
    % Flyer plot, lovingly coded by Paul Pounds, first coded 17/4/02
    % version 2 2004 added scaling and ground display
    % version 3 2010 improved rotor rendering and fixed mirroring bug
    %
    % Displays X-4 flyer position and attitude in a 3D plot.
    % GREEN ROTOR POINTS NORTH
    % BLUE ROTOR POINTS EAST

    % PARAMETERS
    % s defines the plot size in meters
    % swi controls flyer attitude plot; 1 = on, otherwise off.

    % INPUTS
    % 1 Center X position
    % 2 Center Y position
    % 3 Center Z position
    % 4 Yaw angle in rad
    % 5 Pitch angle in rad
    % 6 Roll angle in rad

    % OUTPUTS
    %   None
    ts = [-1 0];

    if ~isfield(vehicle, 'nrotors')
        vehicle.nrotors = 4;    % sensible default for quadrotor function
    end

    switch flag,
        case 0
            [sys,x0,str,ts] = mdlInitializeSizes(ts,plot,enable); % Initialization
        case 3
            sys = mdlOutputs(t,u,s,plot,enable, vehicle); % Calculate outputs
        case {1,2, 4, 9} % Unused flags
            sys = [];
        otherwise
            error(['unhandled flag = ',num2str(flag)]); % Error handling
    end

    % Initialize
function [sys,x0,str,ts] = mdlInitializeSizes(ts,plot,enable)
    % Call simsizes for a sizes structure, fill it in, and convert it
    % to a sizes array.
    sizes = simsizes;
    sizes.NumContStates  = 0;
    sizes.NumDiscStates  = 0;
    sizes.NumOutputs     = 0;
    sizes.NumInputs      = 6;
    sizes.DirFeedthrough = 1;
    sizes.NumSampleTimes = 1;
    sys = simsizes(sizes);
    x0  = [];
    str = [];          % Set str to an empty matrix.
    ts = [0.05 0];

    if enable == 1
        figure(plot);
        clf;
        %colordef(1,'none');
    end
    % End of mdlInitializeSizes.

function sys = mdlOutputs(t,u,s, plot, enable, quad)
    global a1s b1s

    % not quite sure what this is about -- PIC
    if numel(a1s) == [0];
        a1s = zeros(1, quad.nrotors);
        b1s = zeros(1, quad.nrotors);
    end

    % vehicle dimensons
    d = quad.d; %Hub displacement from COG
    r = quad.r; %Rotor radius

    for i = 1:quad.nrotors
        theta = (i-1)/quad.nrotors*2*pi;
        %   Di      Rotor hub displacements (1x3)
        % first rotor is on the x-axis, clockwise order looking down from above
        D(:,i) = [ d*cos(theta); d*sin(theta); 0];
        scal = s(1)/4;
        %Attitude center displacements
        C(:,i) = [ scal*cos(theta); scal*sin(theta); 0];
    end

    if enable == 1
        %draw ground
        figure(plot);
        clf;
        if length(s) == 1
            axis([-s s -s s 0 s]);
        else
            axis([-s(1) s(1) -s(1) s(1) 0 s(2)])
            s = s(1);
        end
        hold on;

        % plot the ground boundaries and the big cross
        plot3([-s -s],[s -s],[0 0],'-b')
        plot3([-s s],[s s],[0 0],'-b')
        plot3([s -s],[-s -s],[0 0],'-b')
        plot3([s s],[s -s],[0 0],'-b')
        plot3([s -s],[-s s],[0 0],'-b')
        plot3([-s s],[-s s],[0 0],'-b')

        %READ STATE
        z = [u(1);u(2);u(3)];
        n = [u(4);u(5);u(6)];

        %PREPROCESS ROTATION MATRIX
        phi = n(1);    %Euler angles
        the = n(2);
        psi = n(3);

        R = [cos(the)*cos(phi) sin(psi)*sin(the)*cos(phi)-cos(psi)*sin(phi) cos(psi)*sin(the)*cos(phi)+sin(psi)*sin(phi);   %BBF > Inertial rotation matrix
            cos(the)*sin(phi) sin(psi)*sin(the)*sin(phi)+cos(psi)*cos(phi) cos(psi)*sin(the)*sin(phi)-sin(psi)*cos(phi);
            -sin(the)         sin(psi)*cos(the)                            cos(psi)*cos(the)];

        %Manual Construction
        %Q3 = [cos(psi) -sin(psi) 0;sin(psi) cos(psi) 0;0 0 1];   %Rotation mappings
        %Q2 = [cos(the) 0 sin(the);0 1 0;-sin(the) 0 cos(the)];
        %Q1 = [1 0 0;0 cos(phi) -sin(phi);0 sin(phi) cos(phi)];
        %R = Q3*Q2*Q1;    %Rotation matrix

        %CALCULATE FLYER TIP POSITONS USING COORDINATE FRAME ROTATION
        F = [1 0 0;0 -1 0;0 0 -1];

        %Draw flyer rotors
        t = [0:pi/8:2*pi];
        for j = 1:length(t)
            circle(:,j) = [r*sin(t(j));r*cos(t(j));0];
        end

        for i = 1:quad.nrotors
            hub(:,i) = F*(z + R*D(:,i)); %points in the inertial frame

            q = 1; %Flapping angle scaling for output display - makes it easier to see what flapping is occurring
            Rr = [cos(q*a1s(i))  sin(q*b1s(i))*sin(q*a1s(i)) cos(q*b1s(i))*sin(q*a1s(i));   %Rotor > Plot frame
                0              cos(q*b1s(i))               -sin(q*b1s(i));
                -sin(q*a1s(i)) sin(q*b1s(i))*cos(q*a1s(i)) cos(q*b1s(i))*cos(q*a1s(i))];

            tippath(:,:,i) = F*R*Rr*circle;
            plot3([hub(1,i)+tippath(1,:,i)],[hub(2,i)+tippath(2,:,i)],[hub(3,i)+tippath(3,:,i)],'b-')
        end

        %Draw flyer
        hub0 = F*z;  % centre of vehicle
        for i = 1:quad.nrotors
            % line from hub to centre plot3([hub(1,N) hub(1,S)],[hub(2,N) hub(2,S)],[hub(3,N) hub(3,S)],'-b')
            plot3([hub(1,i) hub0(1)],[hub(2,i) hub0(2)],[hub(3,i) hub0(3)],'-b')

            % plot a circle at the hub itself
            plot3([hub(1,i)],[hub(2,i)],[hub(3,i)],'o')
        end

        % plot the vehicle's centroid on the ground plane
        plot3([z(1) 0],[-z(2) 0],[0 0],'--k')
        plot3([z(1)],[-z(2)],[0],'xk')

        % label the axes
        xlabel('x');
        ylabel('y');
        zlabel('z (height above ground)');
    end

    sys = [];
    % End of mdlOutputs.

author: Allopart

repo: rbpf-gmapping


text:
function [success] = resetMiji(varargin)
	% This clears Miji from Java's dynamic path and then re-initializes. Use if Miji is not loading normally.
	% Biafra Ahanonu
	% started: 2019.01.28 [14:34:14]
	% inputs
		%
	% outputs
		%
	% usage
		% % Run the following commands one at a time in the command window.
		% resetMiji
		% % An instance of Miji should appear.
		% currP=pwd;Miji;cd(currP);
		% MIJ.exit

	% changelog
		% 2019.09.09 [18:04:10] - Updated to add fallback to modelAddOutsideDependencies to ask for path to Fiji/Miji in case it is not properly in the path or has been moved.
		% 2019.10.15 [12:29:30] - Added flag to prevent recursive loop between resetMiji and modelAddOutsideDependencies.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
        % 2021.12.19 [10:34:00] - Update to use manageMiji to setup Miji, avoid certain issues. Also search for imagej instead of Fiji path.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% DESCRIPTION
	options.nRounds = 1;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	% try
	success = 0;

	for i = 1:options.nRounds
		try
			clear MIJ miji Miji mij;
			javaDyna = javaclasspath('-dynamic');
			% matchIdx = ~cellfun(@isempty,regexpi(javaDyna,'Fiji'));
            matchIdx = ~cellfun(@isempty,regexpi(javaDyna,'imagej'));
			% cellfun(@(x) javarmpath(x),javaDyna(matchIdx));
			javaDynaPathStr = join(javaDyna(matchIdx),''',''');
			if ~isempty(javaDynaPathStr)
				eval(sprintf('javarmpath(''%s'');',javaDynaPathStr{1}))
			else
				disp('Empty Java dynamic path!')
			end
			clear MIJ miji Miji mij;
			% Load Miji so paths added to javaclasspath('-dynamic')
			currP = pwd;
			% Miji;
			% Miji(false);
            ciapkg.io.manageMiji('startStop','setupImageJ');
            MIJ.start;
			cd(currP);
			MIJ.exit;
			% pause(1);
			% java.lang.Runtime.getRuntime().gc;
			% Miji;
			% MIJ.exit;
		catch err
			disp(repmat('@',1,7))
			disp(getReport(err,'extended','hyperlinks','on'));
			disp(repmat('@',1,7))

			try
				modelAddOutsideDependencies('miji','recursionExit',1);
			catch err
				disp(repmat('@',1,7))
				disp(getReport(err,'extended','hyperlinks','on'));
				disp(repmat('@',1,7))
			end
		end
	end

	success = 1;
	% catch err
	% 	success = 0;
	% 	display(repmat('@',1,7))
	% 	disp(getReport(err,'extended','hyperlinks','on'));
	% 	display(repmat('@',1,7))
	% end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function vw = drawROIsMontage(vw)
%
% vw = drawROIsMontage(vw)
%
% Draw the ROIs on a montage view (e.g. several
% flat levels / inplane slices at once).
%
% This function should be called only if
% vw.showROIs is non-zero.  If showROIs=1, draw only the
% selected ROI.  If showROIs=2, draw all the ROIs.
% 09/04 ras, off drawROIs.

% clear line handles for all ROIs (ras 06/06, see selectROI)
for r=1:length(vw.ROIs), vw.ROIs(r).lineHandles = []; end

% set line width parameter
w = 0.5;

s = vw.selectedROI;
n = length(vw.ROIs);

%%%%% get order of ROIs to display
if s
    order = [1:s-1,s+1:n,s];
else
    order = 1:n;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop: draw ROIs according to prefs  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for r = order   % loop through ROIs that are selected for display
    if (r==vw.selectedROI)
        % Selected ROI: now a settable param, ras 05/05
        coords = vw.ROIs(r).coords;
        color = viewGet(vw,'selRoiColor');
    else
        % Non-selected ROI, color from ROIs substruct
        % If showROIs~=2, then set coords=[] so that nothing will be
        % drawn.
        if (vw.ui.showROIs==-2)
            coords = vw.ROIs(r).coords;
        else
            coords = [];
        end
        thisCol=vw.ROIs(r).color;
        % If it's a 'text' color, translate it...
        if (ischar(thisCol))
        nColEntries=length(thisCol);

        for thisColEntry=1:nColEntries
            thisColValue=thisCol(thisColEntry);
            colRGB=[0 0 0];
            switch (thisColValue)
                case 'y', colorRGB = colRGB+[1 1 0];
                case 'm', colorRGB = colRGB+[1 0 1];
                case 'c', colorRGB = colRGB+[0 1 1];
                case 'r', colorRGB = colRGB+[1 0 0];
                case 'g', colorRGB = colRGB+[0 1 0];
                case 'b', colorRGB = colRGB+[0 0 1];
                case 'w', colorRGB = colRGB+[1 1 1];
                case 'k', colorRGB = colRGB+[0 0 0];
                otherwise, colorRGB = [1 1 1];
            end % end switch statement
        end % end loop
            color=colorRGB./nColEntries;
        else
            color=thisCol;
        end

    end

    % Draw the lines around each pixel, w=1/2 because you don't
    % want to connect the centers of the pixels, rather you want to
    % draw around each pixel, e.g, from (x-.5,y-.5) to (x+.5,y-.5).
    if ~isempty(coords)

        % The FLAT view can have a 'rotateImageDegrees' field that
        % specifies a rotation angle for L and R.
        % If this is set, then we have to transform the ROIs by this amount as
        % well to make them register with the anatomy and functional data
        if (strcmp(vw.viewType,'Flat'))
            coords=rotateCoords(vw,coords);
        end

        % convert to coordinates of the montage image
        coords = coords2Montage(vw,coords);
        x = coords(1,:);
        y = coords(2,:);

        % draw boxes around the voxels
        hold on

        % if the view is set to flip the image, then we need to
        % flip the ROI too
        if viewGet(vw, 'flip up down'),
            height = size(vw.ui.image, 1);
            x = height - x + 1;
        end

        for i=1:size(coords,2);
            l(i) = line([y(i)-w, y(i)+w, y(i)+w, y(i)-w, y(i)-w],...
                        [x(i)-w, x(i)-w, x(i)+w, x(i)+w, x(i)-w], ...
                        'Color',color);
        end
        hold off

%         % remember line handles for quick refreshing
%         vw.ROIs(r).lineHandles = l;
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function AnalyzePersistence(Data,AggregationPersistence,LagsSamplAutCorr,Name)

for n=1:length(AggregationPersistence)
    Series_Changes(n).Lag=AggregationPersistence(n);
    SamplAutoCorrs=[];

    for s=1:AggregationPersistence(n)

        Sparse=Data(s:AggregationPersistence(n):end);
        Changes=diff(Sparse);
        SamplAutoCorrs = autocorr(Changes, LagsSamplAutCorr);
        SamplAutoCorrs(1)=[];
    end
    Series_Changes(n).SampleAutocorrelation=mean(SamplAutoCorrs,2);
end

SampleAutocorrelations=[];
for n=1:length(AggregationPersistence)
    SampleAutocorrelations=[SampleAutocorrelations Series_Changes(n).SampleAutocorrelation];
end

figure
bar3(SampleAutocorrelations,'detached')
set(gca,'xlim',[AggregationPersistence(1) AggregationPersistence(end)],'ylim',[1 LagsSamplAutCorr],...
    'zlim',[-.2 .4])
xlabel('aggregation size (days)')
ylabel('lag')
zlabel('autocorrelation')
colormap(.5*(1+gray));
set(gcf,'Name',['        persistence properties of ' Name ' par swap rate'])


text:
function [varargout]=figuremax(varargin)

% function [hf]=figuremax(Cbg,Cdef)
% ------------------------------------------------------------------------
%
% This function opens a maximized figure window using the color settings
% specified.
%
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 30/09/2014
%------------------------------------------------------------------------

%%

switch nargin
    case 0
        Cbg='w';
        Cdef='white';
    case 1
        Cbg=varargin{1};
        Cdef='white';
    case 2
        Cbg=varargin{1};
        Cdef=varargin{2};
    case 3
        Cbg=varargin{1};
        Cdef=varargin{2};
end

figStruct.Color=Cbg; %Figure background color
figStruct.ColorDef=Cdef; %Setting colordefinitions

%Open custom figure
hf=cFigure(figStruct);

%Maximize figure window
maxFig(hf);

if nargout>0
    varargout{1}=hf;
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [fname] = globedems(latlim,lonlim)

    %GLOBEDEMS (30-arc-sec resolution) DEM file names
    %
    % fname = GLOBEDEMS(LATLIM,LONLIM) returns a cellarray of the file
    % names covering the geographic region for GLOBEDEM digital elevation maps.
    % The region is specified by scalar latitude and longitude points, or two
    % element vectors of latitude and longitude limits in units of degrees.
    %
    % The data and some documentation is available over the World-Wide-Web from
    % <http://www.ngdc.noaa.gov/seg/topo/globe.shtml>. The data are available
    % by anonymous ftp from <ftp://ftp.ngdc.noaa.gov/GLOBE_DEM/data/elev/>. The
    % web site also sells copies of the data on CD-ROM.
    %
    % See also: GLOBEDEM

    %  Copyright 1996-2000 Systems Planning and Analysis, Inc. and The MathWorks, Inc.
    %  Written by:  A. Kim, W. Stumpf, L. Job
    %  $Revision: 1399 $    $Date: 2006-08-11 11:19:27 +0200 (Fr, 11 Aug 2006) $

    if nargin~=2
        error('Incorrect number of arguments')
    end

    % ensure row vectors
    latlim = latlim(:)';
    lonlim = lonlim(:)';

    if  isequal(size(latlim),[1 1])
        latlim = latlim*[1 1];
    elseif ~isequal(size(latlim),[1 2])
        error('Latitude limit input must be a scalar or 2 element vector')
    end

    if isequal(sort(size(lonlim)),[1 1])
        lonlim = lonlim*[1 1];
    elseif ~isequal(sort(size(lonlim)),[1 2])
        error('Longitude limit input must be a scalar or 2 element vector')
    end

    % read names and bounding rectangle limits

    [fnames,YMIN,YMAX,XMIN,XMAX,rtile,ctile] = textread('globedems.dat','%s %d %d %d %d %d %d');

    % case where dateline is not crossed
    if lonlim(1) <= lonlim(2)
        do = ...
            find( ...
            (...
            (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
            (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
            (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
            (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
            ) ...
            &...
            (...
            (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
            (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
            (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
            (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
            )...
            );
    end

    % case where the dateline is crossed
    if lonlim(1) > lonlim(2)
        lmin = lonlim(1); lmax = lonlim(2);
        lonlim(2) = 180;
        % do eastern side of the dateline first
        doEAST = ...
            find( ...
            (...
            (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
            (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
            (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
            (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
            ) ...
            &...
            (...
            (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
            (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
            (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
            (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
            )...
            );
        % do western side of the dateline second
        lonlim(1) = -180; lonlim(2) = lmax;
        doWEST = ...
            find( ...
            (...
            (latlim(1) <= YMIN & latlim(2) >= YMAX) | ... % tile is completely within region
            (latlim(1) >= YMIN & latlim(2) <= YMAX) | ... % region is completely within tile
            (latlim(1) >  YMIN & latlim(1) <  YMAX) | ... % min of region is on tile
            (latlim(2) >  YMIN & latlim(2) <  YMAX)   ... % max of region is on tile
            ) ...
            &...
            (...
            (lonlim(1) <= XMIN & lonlim(2) >= XMAX) | ... % tile is completely within region
            (lonlim(1) >= XMIN & lonlim(2) <= XMAX) | ... % region is completely within tile
            (lonlim(1) >  XMIN & lonlim(1) <  XMAX) | ... % min of region is on tile
            (lonlim(2) >  XMIN & lonlim(2) <  XMAX)   ... % max of region is on tile
            )...
            );
        % concatenate indices
        do = [doEAST doWEST];
    end

    if ~isempty(do)
        fname = fnames(do);
    else
        fname = [];
    end

author: CelsoReyes


text:
%DEMO_REGRESSION_SPARSE2  Regression demo comparing different sparse
%                         approximations with optimization of inducing
%                         variables
%
%  Description
%    A regression problem with one input variable and one output
%    variable with Gaussian noise. The output is assumed to be
%    realization of additive functions and Gaussian noise.
%
%    For standard full GP demonstration, see for example
%    DEMO_REGRESSION1, and for detailed discussion, Rasmussen and
%    Williams (2006). For more sparse demonstrations including use
%    of compact support covariance functions see DEMO_REGRESSION2,
%    DEMO_MODELASSESMENT2, and DEMO_SPARSEAPPROX.
%
%    In this demo, sparse approximations for the full GP model are
%    compared. We use
%      - FIC, fully independent conditional
%      - DTC, deterministic training conditional
%      - VAR, variational approach
%    For illustration purposes the hyperparameters from the full GP
%    are used for the sparse models and only the inducing variables
%    are optimised.
%
%    For technical details, see Quinonero-Candela and Rasmussen
%    (2005) for the FIC and DTC models and Titsias (2009) for the
%    VAR model.
%
%    We use a simple one dimensional data set to present the three
%    methods.
%
%  See also DEMO_REGRESSION1, DEMO_REGRESSION2, DEMO_REGRESSION_SPARSE1
%
%
%  References:
%
%    Quinonero-Candela, J. and Rasmussen, C. E. (2005). A Unifying
%    View of Sparse Approximate Gaussian Process Regression. Journal
%    of Machine Learning Research.
%
%    Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian
%    Processes for Machine Learning. The MIT Press.
%
%    Titsias, M. K. (2009). Variational Model Selection for Sparse
%    Gaussian Process Regression. Technical Report, University of
%    Manchester.

% Copyright (c) 2010 Heikki Peura, 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.

% Set randomstream for reproducing same results
prevstream=setrandstream();

% Start by creating 1D data
xx=linspace(1,10,901);

% Choose a subset of data so that the data are less dense in the right end.
% xt are the inputs and yt are the outputs, xstar are the values we want to
% predict.
x1=logspace(0,1,100);
x1=round(x1*100)-99;
x=xx(x1)';
y=2*sin(4*x)+0.2*randn(size(x));
xt=[1:0.01:14]';
[n,nin] = size(x);

fprintf('Full GP\n')
% Initialize full GP with a squared exponential component and set
% priors for their parameters.
pl = prior_t('s2', 1);
pm = prior_logunif();
pn = prior_logunif();

gpcfse = gpcf_sexp('lengthScale',0.5,'magnSigma2',1, 'lengthScale_prior', pl, 'magnSigma2_prior', pm);
lik = lik_gaussian('sigma2', 0.1, 'sigma2_prior', pn);

gp = gp_set('lik', lik, 'cf', gpcfse, 'jitterSigma2', 1e-6);

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

[Eft_full, Varft_full] = gp_pred(gp, x, y, xt);
Varft_full = Varft_full + gp.lik.sigma2;

figure
% Blue crosses are the initial inducing input locations, red ones are
% the optimised ones. Black circles represent the distance to the next
% optimized location, with a dashed trendline.');

subplot(2,2,1);hold on;
plot(xt,Eft_full,'k', 'LineWidth', 2)
plot(xt,Eft_full-2.*sqrt(Varft_full),'--','Color',[0 0.5 0])
plot(xt,Eft_full+2.*sqrt(Varft_full),'--','Color',[0 0.5 0])
plot(x,y,'.', 'MarkerSize',7)
ylim([-3 3])
title('FULL GP')

fprintf('FIC GP\n')
% Run FIC approximation for the same data: choose the inducing
% inputs Xu, then proceed with the inference with the optimized
% parameters from the full GP: here, we optimize only the locations
% of the inducing inputs for the FIC model.
Xu=round(10+90*rand(18,1))/10; % Random placement

% Change type to FIC, add inducing inputs, and optimize only inducing inputs
gp_fic = gp_set(gp, 'type','FIC','X_u',Xu,'infer_params','inducing');

% Set the options for the optimization
opt=optimset('TolFun',1e-4,'TolX',1e-4);
% Optimize with the scaled conjugate gradient method
gp_fic=gp_optim(gp_fic,x,y,'opt',opt);

[Eft_fic, Varft_fic] = gp_pred(gp_fic, x, y, xt);
Varft_fic = Varft_fic + gp_fic.lik.sigma2;

subplot(2,2,2);hold on;
h1=plot(xt,Eft_fic,'k', 'LineWidth', 2);
h2=plot(xt,Eft_fic-2.*sqrt(Varft_fic),'--','Color',[0 0.5 0]);
plot(xt,Eft_fic+2.*sqrt(Varft_fic),'--','Color',[0 0.5 0])
h3=plot(x,y,'.', 'MarkerSize',7);
h4=plot(Xu, -2.8, 'bx', 'MarkerSize', 5, 'LineWidth', 2);
h5=plot(gp_fic.X_u, -3, 'rx', 'MarkerSize', 5, 'LineWidth', 2);
legend([h1  h2 h3 h4(1) h5(1)],'Ef','95% CI','Data','Initial X_u','Optimized X_u')
% plot diff of sorted X_u and regress line for that
%XuSorted=sort(gp_fic.X_u);
%dXuSorted=diff(XuSorted);
%bb=regress(dXuSorted,[ones(size(dXuSorted)) XuSorted(1:end-1)]);
%plotbb=bb(1)+(min(XuSorted):0.1:max(XuSorted))*bb(2);
%plot(XuSorted(1:end-1),dXuSorted,'ko');
%plot(min(XuSorted):0.1:max(XuSorted),plotbb,'k--')
ylim([-3 3])
title('FIC')

fprintf('VAR GP\n')
% Run the VAR model similarly to the FIC model with the same
% starting inducing inputs. The difference in the optimized results
% is notable. The VAR model places the inducing inputs quite evenly
% (slightly increasing as the data becomes more sparse), with
% predictions closely matching the full GP model. The other two
% sparse approximations yield less reliable results.
gp_var = gp_set(gp,'type','VAR','X_u',Xu,'infer_params','inducing');

% Set the options for the optimization
opt=optimset('TolFun',1e-4,'TolX',1e-4);
% Optimize with the scaled conjugate gradient method
gp_var=gp_optim(gp_var,x,y,'opt',opt);

[Eft_var, Varft_var] = gp_pred(gp_var, x, y, xt);
Varft_var = Varft_var + gp_var.lik.sigma2;

subplot(2,2,4);hold on
plot(xt,Eft_var,'k', 'LineWidth', 2)
plot(xt,Eft_var-2.*sqrt(Varft_var),'--','Color',[0 0.5 0])
plot(xt,Eft_var+2.*sqrt(Varft_var),'--','Color',[0 0.5 0])
plot(x,y,'.', 'MarkerSize',7)
plot(gp_var.X_u, -3, 'rx', 'MarkerSize', 5, 'LineWidth', 2)
plot(Xu, -2.8, 'bx', 'MarkerSize', 5, 'LineWidth', 2)
% plot diff of sorted X_u and regress line for that
%XuSorted=sort(gp_var.X_u);
%dXuSorted=diff(XuSorted);
%bb=regress(dXuSorted,[ones(size(dXuSorted)) XuSorted(1:end-1)]);
%plotbb=bb(1)+(min(XuSorted):0.1:max(XuSorted))*bb(2);
%plot(XuSorted(1:end-1),dXuSorted,'ko');
%plot(min(XuSorted):0.1:max(XuSorted),plotbb,'k--')
ylim([-3 3])
title('VAR')

fprintf('DTC GP\n')
% Run the DTC model similarly to the FIC model with the same starting
% inducing inputs. The difference in the optimized results is notable.
gp_dtc = gp_set(gp,'type','DTC','X_u',Xu,'infer_params','inducing');

% Set the options for the optimization
opt=optimset('TolFun',1e-4,'TolX',1e-4);
% Optimize with the scaled conjugate gradient method
gp_dtc=gp_optim(gp_dtc,x,y,'opt',opt);

[Eft_dtc, Varft_dtc] = gp_pred(gp_dtc, x, y, xt);
Varft_dtc = Varft_dtc + gp_dtc.lik.sigma2;

subplot(2,2,3);hold on
plot(xt,Eft_dtc,'k', 'LineWidth', 2)
plot(xt,Eft_dtc-2.*sqrt(Varft_dtc),'--','Color',[0 0.5 0])
plot(xt,Eft_dtc+2.*sqrt(Varft_dtc),'--','Color',[0 0.5 0])
plot(x,y,'.', 'MarkerSize',7)
plot(gp_dtc.X_u, -3, 'rx', 'MarkerSize', 5, 'LineWidth', 2)
plot(Xu, -2.8, 'bx', 'MarkerSize', 5, 'LineWidth', 2)
% plot diff of sorted X_u and regress line for that
%XuSorted=sort(gp_dtc.X_u);
%dXuSorted=diff(XuSorted);
%bb=regress(dXuSorted,[ones(size(dXuSorted)) XuSorted(1:end-1)]);
%plotbb=bb(1)+(min(XuSorted):0.1:max(XuSorted))*bb(2);
%plot(XuSorted(1:end-1),dXuSorted,'ko');
%plot(min(XuSorted):0.1:max(XuSorted),plotbb,'k--')
ylim([-3 3])
title('DTC')

% Set back initial random stream
setrandstream([],prevstream);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function symmetric_sparse_size_test ( )

%*****************************************************************************80
%
%% SYMMETRIC_SPARSE_SIZE_TEST tests SYMMETRIC_SPARSE_SIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 September 2012
%
%  Author:
%
%    John Burkardt.
%
%  Local parameters:
%
%    Local, integer D, the spatial dimension.
%
%    Local, integer MAXK, the maximum level to check.
%
  test_num = 3;

  dim_test = [ 5, 5, 3 ];
  nodes1 = [ ...
   0.0, 0.0, 0.0, 0.0, 0.0, 1.0;
   0.0, 0.0, 0.0, 0.0, 1.0, 0.0;
   0.0, 0.0, 0.0, 1.0, 0.0, 0.0;
   0.0, 0.0, 1.0, 0.0, 0.0, 0.0;
   0.0, 1.0, 0.0, 0.0, 0.0, 0.0 ]';
  nodes2 = [ ...
    0.0, 0.0, 0.0, 0.0, 0.0, ...
    0.0, 0.0, 0.0, 0.0, 0.0, ...
    0.0, 0.0, 0.0, 0.0, 0.0, ...
    1.0, 1.0, 1.0, 1.0, 1.0, ...
    1.73205; ...
    0.0, 0.0, 0.0, 0.0, 0.0, ...
    0.0, 0.0, 0.0, 0.0, 0.0, ...
    1.0, 1.0, 1.0, 1.0, 1.73205, ...
    0.0, 0.0, 0.0, 0.0, 1.0, ...
    0.0; ...
    0.0, 0.0, 0.0, 0.0, 0.0, ...
    0.0, 1.0, 1.0, 1.0, 1.73205, ...
    0.0, 0.0, 0.0, 1.0, 0.0, ...
    0.0, 0.0, 0.0, 1.0, 0.0, ...
    0.0; ...
    0.0, 0.0, 0.0, 1.0, 1.0, ...
    1.73205, 0.0, 0.0, 1.0, 0.0, ...
    0.0, 0.0, 1.0, 0.0, 0.0, ...
    0.0, 0.0, 1.0, 0.0, 0.0, ...
    0.0; ...
    0.0, 1.0, 1.73205, 0.0, 1.0, ...
    0.0, 0.0, 1.0, 0.0, 0.0, ...
    0.0, 1.0, 0.0, 0.0, 0.0, ...
    0.0, 1.0, 0.0, 0.0, 0.0, ...
    0.0 ]';
  nodes3 = [ ...
    0.0, 0.0, 0.0, 0.0, 0.0, ...
    0.0, 0.0, 0.0, 0.0, 0.0, ...
    0.0, 0.0, 0.741964, 1.0, 1.0, ...
    1.0, 1.0, 1.0, 1.0, 1.73205, ...
    1.73205, 1.73205, 2.33441; ...
    0.0, 0.0, 0.0, 0.0, 0.0, ...
    0.741964, 1.0, 1.0, 1.0, 1.73205, ...
    1.73205, 2.33441, 0.0, 0.0, 0.0, ...
    0.0, 1.0, 1.0, 1.73205, 0.0, ...
    0.0, 1.0, 0.0; ...
    0.0, 0.741964, 1.0, 1.73205, 2.33441, ...
    0.0, 0.0, 1.0, 1.73205, 0.0, ...
    1.0, 0.0, 0.0, 0.0, 1.0, ...
    1.73205, 0.0, 1.0, 0.0, 0.0, ...
    1.0, 0.0, 0.0 ]';
  r_test = [ 6, 21, 23 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SYMMETRIC_SPARSE_SIZE_TEST\n' );
  fprintf ( 1, '  Given a symmetric sparse grid rule represented only by\n' );
  fprintf ( 1, '  the points with positive values, determine the total number\n' );
  fprintf ( 1, '  of points in the grid.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For dimension DIM, we report\n' );
  fprintf ( 1, '  R, the number of points in the positive orthant, and\n' );
  fprintf ( 1, '  R2, the total number of points.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       DIM         R        R2\n' );
  fprintf ( 1, '\n' );

  x0 = 0.0;

  for test = 1 : test_num

    r = r_test(test);
    dim = dim_test(test);

    if ( test == 1 )
      r2 = symmetric_sparse_size ( r, dim, nodes1, x0 );
    elseif ( test == 2 )
      r2 = symmetric_sparse_size ( r, dim, nodes2, x0 );
    elseif ( test == 3 )
      r2 = symmetric_sparse_size ( r, dim, nodes3, x0 );
    end

    fprintf ( 1, '  %8d  %8d  %8d\n', dim, r, r2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [lfp, spike, stm, bhv] = spass2fieldtrip(dirname, varargin)

% SPASS2FIELDTRIP reads data from a set of SPASS data files and converts
% the contents into data structures that FieldTrip understands. Note that
% dependent on the SPASS data it might be required to change some
% hard-coded parameters inside this function.
%
% Use as
%   [lfp, spike, stm, bhv] = spass2fieldtrip(dirname)
% Optionally you can specify the sample rate as key-value pairs
%  'fsample_ana' - default 1000
%  'fsample_swa' - default 32000
%
% The specified directory should contain the SPASS files, and the files should have
% the same name as the directory.
%
% The swa and sti input file are combined into the spike output structure.
% For the rest of the data it is trivial how the input and output relate.
%
% For example, if you specify
%   [lfp, spike, bhv, stm] = spass2fieldtrip('jeb012a02')
% then the following files should exist:
%   'jeb012a02/jeb012a02.ana'
%   'jeb012a02/jeb012a02.swa'
%   'jeb012a02/jeb012a02.spi'
%   'jeb012a02/jeb012a02.stm'
%   'jeb012a02/jeb012a02.bhv'
%
% Subsequently you can analyze the data in FieldTrip, or write the spike
% waveforms to a nex file for offline sorting using
%   ft_write_spike('jeb012a02_ch1.nex', spike, 'dataformat', 'plexon_nex', 'chanindx', 1)
%   ft_write_spike('jeb012a02_ch2.nex', spike, 'dataformat', 'plexon_nex', 'chanindx', 2)
%   ft_write_spike('jeb012a02_ch3.nex', spike, 'dataformat', 'plexon_nex', 'chanindx', 3)
%
% See also NUTMEG2FIELDTRIP, LORETA2FIELDTRIP, FIELDTRIP2SPSS

% Copyright (C) 2007, Robert Oostenveld
%
% $Id$

fsample_ana = ft_getopt(varargin, 'fsample_ana', 1000);
fsample_swa = ft_getopt(varargin, 'fsample_swa', 32000);

anafile = fullfile(dirname, [dirname '.ana']);
swafile = fullfile(dirname, [dirname '.swa']);
spifile = fullfile(dirname, [dirname '.spi']);
stmfile = fullfile(dirname, [dirname '.stm']);
bhvfile = fullfile(dirname, [dirname '.bhv']);

if ~exist(anafile, 'file'), ft_error('the file "%s" does not exist', anafile); end
if ~exist(swafile, 'file'), ft_error('the file "%s" does not exist', swafile); end
if ~exist(spifile, 'file'), ft_error('the file "%s" does not exist', spifile); end
if ~exist(stmfile, 'file'), ft_error('the file "%s" does not exist', stmfile); end
if ~exist(bhvfile, 'file'), ft_error('the file "%s" does not exist', bhvfile); end

% read the data
fprintf('reading %s\n', anafile); ana = read_labview_dtlg(anafile, 'int16');
fprintf('reading %s\n', swafile); swa = read_labview_dtlg(swafile, 'int16');
fprintf('reading %s\n', spifile); spi = read_labview_dtlg(spifile, 'int32');
fprintf('reading %s\n', stmfile); stm = read_labview_dtlg(stmfile, 'int32');
fprintf('reading %s\n', bhvfile); bhv = read_labview_dtlg(bhvfile, 'int32');

% determine the number of trials
ntrials = numel(bhv.data{1});

nchans = numel(ana.data)./ntrials;
ana.data = reshape(ana.data, [nchans, ntrials]);

% prepare the continuous LFP data
lfp         = [];
lfp.trial   = {};
lfp.time    = {};
lfp.label   = {};
lfp.fsample = fsample_ana;
for i=1:ntrials
  tmp = cell2mat(ana.data(:,i)')';
  lfp.trial{i} = tmp;
  lfp.time{i}  = ((1:size(tmp,2))-1)/fsample_ana;
  nsamples(i)  = length(lfp.time{i});
end
for i=1:size(ana.data,1)
  lfp.label{i,1} = sprintf('chan%d', i);
end

% the data is trial based, try to estimate the time between subsequent
% trials, or better: the time between subsequent stimuli
isi = 10^ceil(log10(max(nsamples)+1));

nchans = numel(swa.data)./ntrials;
swa.data = reshape(swa.data, [nchans, ntrials]);
spi.data = reshape(spi.data, [nchans, ntrials]);

% prepare the spike data
spike           = [];
spike.label     = {};
spike.waveform  = {}; % 1xnchans cell-array, each element contains a matrix (Nsamples X Nspikes), can be empty
spike.timestamp = {}; % 1xnchans cell-array, each element contains a vector (1 X Nspikes)
spike.unit      = {}; % 1xnchans cell-array, each element contains a vector (1 X Nspikes)
for i=1:size(swa.data,1)
  spike.label{i} = sprintf('chan%d', i);
end

nspikes = zeros(nchans, ntrials);
for i=1:nchans
  for j=1:ntrials
    nspikes(i,j) = length(swa.data{i,j});
  end
end

% reinsert the inter-trial intervals, this links both the LFP and the spike
% timestamps to a common continuous timeaxis
for j=1:nchans
  for i=1:ntrials
    spi.data{j,i} = spi.data{j,i} + (i-1)*(fsample_swa./fsample_ana)*isi;
  end
end

% convert the spike timestamps and waveforms to a FieldTrip-compatible format
for i=1:nchans
  spike.waveform{i}  = cell2mat(swa.data(i,:));
  spike.timestamp{i} = cell2mat(spi.data(i,:)')';
  spike.unit{i}      = ones(size(spike.timestamp{i}));  % since unsorted
end

% prepare the other output
stm = stm.data{1}(:);
bhv = bhv.data{1}(:);

% remember where the lfp trials are on the imaginary continuous timeaxis,
% this links both the LFP and the spike timestamps to a common continuous
% timeaxis
trl = zeros(ntrials, 3);
for i=1:ntrials
  begsample = (i-1)*isi+1;
  endsample = begsample+nsamples(i)-1;
  offset    = 0;
  trl(i,:) = [begsample endsample offset];
end

% store the header information
lfp.hdr.FirstTimeStamp = 0;
lfp.hdr.TimeStampPerSample = fsample_swa./fsample_ana;

author: spm


text:
function [ u_i ] = rotate_b2i( u_b, phi, theta, psi )
% ROTATE_B2I - Rotation of an object (vector/matrix) from the body to the
% inertial 3D-space
%
%
%

    cr = cos(phi);
    cp = cos(theta);
    cy = cos(psi);
    sr = sin(phi);
    sp = sin(theta);
    sy = sin(psi);

    % Rotation matrix from inertial frame to body frame

    Rbi = [cp*cy, sr*sp*cy-cr*sy, cr*sp*cy+sr*sy;...
            cp*sy, sr*sp*sy+cr*cy, cr*sp*sy-sr*cy;...
            -sp,   sr*cp,          cr*cp];

    u_i = Rbi*u_b;

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function M = oned_bilinear ( kernel, phi, test, w_g )

%*****************************************************************************80
%
%% ONE_BILINEAR evaluates a bilinear form.
%
%  Discussion:
%
%    This function computes integral ( kernel * phi * test }
%
%    The calculation that is carried out is equivalent to:
%
%      [ n_gauss, n_test ] = size ( test );
%      [ n_gauss, n_dof ] = size ( phi );
%
%      M = zeros ( n_test, n_dof );
%      for i = 1 : n_test
%        for j = 1 : n_dof
%          M(i,j) = ( kernel' .* test(:,i)' ) * ( phi(:,j) .* w_g );
%        end
%      end
%
%  Parameters:
%
%    Input, real KERNEL(N_GAUSS), the kernel function in the integral,
%    evaluated at the Gauss points.
%
%    Input, real PHI(N_GAUSS,N_DOF), the element test functions evaluated
%    at the Gauss points.
%
%    Input, real TEST(N_GAUSS,N_TEST), the test functions evaluated at the
%    Gauss points.
%
%    Input, real W_G(N_GAUSS,1), a column vector of Gauss weights.
%
%    Output, real M(N_TEST,N_DOF), the bilinear form.
%
  M = test' * diag ( kernel .* w_g ) * phi;

  return
end


text:
%% Potts-denoising, Gaussian noise

%%
% Load signal
original = loadPcwConst('sample1');
% Add Gaussian noise with parameter sigma = 0.1
f = original + 0.1 * randn(size(original));

%%
% Compute $L^2$ Potts estimator
figure(1)
pottsL2 = minL2Potts(f, 0.1);
showPotts(f, pottsL2, original, 'L2-Potts')

%%
% Compute $L^1$ Potts estimator
figure(2)
pottsL1 = minL1Potts(f, 0.4);
showPotts(f, pottsL1, original, 'L1-Potts')

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function I = integral2(f, S)
%INTEGRAL2    Surface integral of a CHEBFUN3.
%   INTEGRAL2(F, S) returns integral of the CHEBFUN3 object F over the
%   parametric surface S defined as a CHEBFUN2V object.
%
%   I = INTEGRAL2(F) is the same as I = SUM2(F).
%
% See also CHEBFUN3/INTEGRAL, CHEBFUN3/INTEGRAL3, CHEBFUN3/SUM,
% CHEBFUN3/SUM2 and CHEBFUN3/SUM3.

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

% Developer Note: If F = F(x,y,z) is a CHEBFUN3 and
% S = {(u,v)\in DOM, s.t. x = x(u,v), y = y(u,v) and z = z(u,v)}, is a
% parametric surface represented by a CHEBFUN2V object, then
% \int \int_S F(x,y,z) dS = \int \int_DOM F(x(u,v), y(u,v), z(u,v))) ...
%                                                 norm(cross(r_u, r_v)) dA.
% Note that the domain of f should contain the range of S.
% TODO: Check for this in the code.

% Empty check:
if ( isempty(f) )
    I = [];
    return
end

if ( nargin == 1 )
    % Double definite integral:
    I = sum2(f);

   elseif ( nargin == 2 )
       if ( isa(S, 'chebfun2v') )
           % Integral over a parametric surface represented by a CHEBFUN2V.
           % Get the surface:
           S_compon = S.components;
           S1 = S_compon{1};
           S2 = S_compon{2};
           S3 = S_compon{3};

           % Surface integral:
           diffCu = diff(S, 1, 2); % Note the Chebfun2 convention to
           diffCv = diff(S, 1, 1); % use 2 for the 1st variable.
           ds = cross(diffCu, diffCv);
           op = @(u,v) feval(f, feval(S1, u, v), feval(S2, u, v), ...
               feval(S3, u, v));
           I = sum2(chebfun2(op, S1.domain).*ds);

       elseif ( isvector(S) && numel(S) == 2 )
           % Double definite integral over specified dimensions:
           I = sum2(f, S);
       end
else
    error('CHEBFUN:CHEBFUN3:integral2:nargin', ['Incorrect number of '...
        'input arguments.']);
end

end

function ds = cross(F, G)
H = [F(2).*G(3) - F(3).*G(2);
     F(3).*G(1) - F(1).*G(3);
     F(1).*G(2) - F(2).*G(1)];
% Developer note: In principle, here we should use
% ds = sqrt(H(1).^2 + H(2).^2 + H(3).^2);
% which uses chebfun2/sqrt. But, that code calls a singleSingTest
% subroutine which sometimes gives error even in this case where the input
% to sqrt is always nonnegative. We bypass that code by calling chebfun2
% constructor as follows:
ds = chebfun2(@(u,v) sqrt(abs(feval(H(1),u,v).^2 + feval(H(2), u, v).^2 +...
    feval(H(3), u, v).^2)), H(1).domain);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function jdoric( n,mink,maxk,a,b )
%RIC :
%   Generates restricted and unrestricted integer compositions
%   n = integer to be partitioned
%   kmin = min. no. of summands
%   kmax = max. no. of summands
%   a = min. value of summands
%   b = max.value of summands
%
% from : "A Unified Approach to Algorithms Generating Unrestricted
%           and Restricted Integer Compositions and Integer Partitions"
%   J. D. OPDYKE, J. Math. Modelling and Algorithms (2009) V9 N1, p.53 - 97
%
% Matlab implementation :
% Theophanes E. Raptis, DAT-NCSRD 2010
% http://cag.dat.demokritos.gr
% rtheo@dat.demokritos.gr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cell = [];
rowdec = 0;
for i=mink:maxk
    in = n/i;
    if a>1 rowdec = i; end
    if a<=in && in <= b N2N(n,i,a,b,n-1-rowdec,i-1,0,0,cell); end
end
end

function N2N(n,i,a,b,row,col,level,cumsum,cell)
if col~=0
    if col==1
        jmax = max(a,n-cumsum-b);
        jmin = min(b,n-cumsum-a);
        for j=jmax : jmin
            cell(i-1) = j; cell(i) = n-cumsum-j;
            disp(cell);
        end
    else
        cell(level+1) = a;
        tmp = cumsum + a;
        ntmp = round((n - tmp)/(i-level-1));
        if a <= ntmp && ntmp <= b && cell(level+1)
            N2N(n,i,a,b,row-a+1,col-1,level+1,tmp,cell);
        else
            for q=1:min((b-a),(row-a)-(col-1))
                cell(level+1) = cell(level+1)+1;
                tmp = tmp + 1;
                ntmp = round((n - tmp)/(i-level-1));
                if a <= ntmp && ntmp <= b && cell(level+1)
                    q2 = q; q = min((b-a),(row-a)-(col-1));
                    N2N(n,i,a,b,row-a-q2+1,col-1,level+1,tmp,cell);
                end
            end
        end
    end
else disp(n)
end
if level>0 && row>1
    cell(level) = cell(level)+1;
    cumsum = cumsum + 1;
    if cell(level)<a
        cumsum = cumsum - cell(level) + a;
        cell(level) = a;
        row = row + cell(level) - a;
    end
    toploop = min(b-cell(level),row-col-1);
    npart = round((n-cumsum)/(i-level));
    if a<=npart && npart<=b && cell(level)<=b
        N2N(n,i,a,b,row-1,col,level,cumsum,cell);
    else
        for p=1:toploop
            cell(level) = cell(level)+1;
            cumsum = cumsum + 1;
            npart = round((n-cumsum)/(i-level));
            if a<=npart && npart<=b && cell(level)<=b
                p2 = p; p = toploop;
                N2N(n,i,a,b,row-p2,col,level,cumsum,cell);
            end
        end
    end
end
end


text:
function [clusters,subclusters] = cluster_princomp(clusters,varargin)
% function [clusters,subclusters] = cluster_princomp(clusters,[behavioral score vector],[corr flag],[plotflag],[locflag])
%
% ALSO TRY: subcluster_montage(subclusters{1}) % to plot the output
%
% clusters is structure of clusters from tor_extract_rois.m
% behavioral vector is row vector of behavioral or other scores to correlate
% corr flag:  *1 = work on correlations among voxels, 2 = work on covariance
% plotflag:   *1 = yes, 0 = no.  plots.
% locflag:    1 yes, *0 no; add XYZ voxel locations (scaled) to data submitted to clustering
%             pushes voxels closer in space to be classified in the same cluster
%
% try this to test the program on random data:
% cl(1).all_data = randn(23,30);cl(1).numVox = 30;cl = cluster_princomp(cl,EXPT.behavior,1,1);
% cl(1).all_data(:,1:10) = cl(1).all_data(:,1:10) + 10; cl = cluster_princomp(cl,EXPT.behavior,1,1);
% cl(1).all_data(:,25:30) = cl(1).all_data(:,25:30) + repmat((EXPT.behavior .* 3)',1,6);
% cl(1).all_data(:,21:24) = cl(1).all_data(:,21:24) + repmat((1:23)',1,4);
% cl = cluster_princomp(cl,EXPT.behavior,1,1);
% mean-center everything now:
% cl.PCA = []; cl.all_data - cl.all_data - repmat(mean(cl.all_data),size(cl.all_data,1),1);
% cl = cluster_princomp(cl,EXPT.behavior,1,1);
% add another correlated group:
% cl.all_data(:,1:5) = cl.all_data(:,1:5) + repmat(rand(23,1)*5,1,5);
% cl = cluster_princomp(cl,EXPT.behavior,1,1);
%
% if component scores are used and correlated with behavior, this means that the subjects
% tend to show the behavioral effect who also show the pattern associated with comp. x.
% this may mean high on a number of voxels, or high on some and low on others.
% the weights may be used to interpret what the components mean, and this can be done
% graphically.
%
% t-tests on component scores have ambiguous interpretations, because a high t-score
% may indicate negative values or close-to-zero values on some voxels.
% a component could have the interpretation, "high on this component means high on V1
% and low on V2."
%
% classifying voxels is done using cluster analysis (hierarchical, centroid linkage)
% on the voxels (observations) using the PCA weights (eigenvectors) as variables.
% This lets the clustering algorithm work in the reduced variable space with dimensionality
% equal to the number of components.
% The max number of clusters is restricted based on the gradient of the eigenvalues in the PCA
% maxclusters = 1 + the number of eigenvalues with gradient at least 20% of the initial drop
% from 1 to 2 eigenvalues.
%
% Requires clustering library in Matlab.
% Robust option also uses the robust PCA algorithm RAPCA,
% created by:
% Hubert, M., Rousseeuw, P.J., Verboven, S. (2002),
%  "A fast method for robust principal components with applications to chemometrics", by Mia Hubert, Peter J. Rousseeuw,
%  Chemometrics and Intelligent Laboratory Systems, 60, 101-111.
%
%

corrflag = 1; plotflag = 1; robustflag = 1; locflag = 0;
if length(varargin) > 1, corrflag = varargin{2};, end
if length(varargin) > 2, plotflag = varargin{3};, end
if length(varargin) > 3, locflag = varargin{4};, end

for i = 1:length(clusters)
    subclusters{i} = []; clusters(i).PCA = [];

    % check for NaNs and remove those voxels
    tst = any(isnan(clusters(i).all_data),1);
    clusters(i).all_data(:,tst) = [];
    clusters(i).XYZ(:,tst) = [];
    clusters(i).XYZmm(:,tst) = [];
    clusters(i).Z(:,tst) = [];
    if any(tst),disp(['Warning! Removed ' num2str(sum(tst)) ' voxels with NaN values.']),end

    a = clusters(i).all_data;

    % scale here, because robust pca doesn't use correlations, does it?
    % but robust PCA seems to be unaffected by scale changes on some variables
    if corrflag, a = scale(a);,end

    % if we choose to add the XYZ flag to add locations to clustering criteria
    if locflag,
        wfactor = round(size(a,1) ./ 3);                       % weight for loc; higher = more weight on location
        xyztmp = repmat(scale(clusters(i).XYZ')',wfactor,1);   % scale to make comparable to img values,
                                                                % but multiply by weighting factor
        a = [a; xyztmp];
    end

    if size(a,2) > 2    % must have 3 voxels to try clustering

    % -------------------------------------------------------------------------------
    % * All the real work is done here.  Compute pc's
    % -------------------------------------------------------------------------------

    if ~robustflag
        [clusters(i).PCA.pcomps,clusters(i).PCA.weights,clusters(i).PCA.eigval,clusters(i).PCA.class] = pc(a,corrflag);

        % automatically pick number of clusters, based on gradient in eigenvalues
        g = abs(gradient(clusters(i).PCA.eigval));
        maxclusters = sum(g > g(1).*.2) + 1;
    else

        % pick number of dimensions by hand
        disp(' ')
        disp([num2str(clusters(i).numVox) ' voxels in main cluster'])
        fprintf(1,'%3.2f Observations on %3.2f voxels\n',size(a,1),size(a,2))
        fprintf(1,'Save at least 2 eigenvectors to do clustering');
        % doing this on a' means voxels are observations, conditions/subj scores are variables
        % we'll use the scores, which has voxels as rows and components as columns, to classify
        out = rapca(a');
        clusters(i).PCA.pcomps = out.T;
        clusters(i).PCA.weights = out.P;
        clusters(i).PCA.eigval = out.L;
        maxclusters = length(out.L);
    end

    % -------------------------------------------------------------------------------
    % * All the real work is done here.  Classify
    % -------------------------------------------------------------------------------

    if size(clusters(1).all_data,1) > 12,
        disp('More than 12 dimensions (observations per voxel) in original data - using eigenvectors to classify')
        close all; try, pack, catch, end
        % we pick the number of CLASSES separately by hand, because # components not a good indicator
        % of how many classes there are
        clusters(i).PCA.class = docluster(out.T',[],plotflag);
    else
        clusters(i).PCA.class = docluster(out.T',[],plotflag);
    end

    % clean up if using locflag
    clusters(i).PCA.locflag = locflag;
    if locflag,
        clusters(i).PCA.pcomps = clusters(i).PCA.pcomps(1:size(clusters(i).all_data,1),:);
        %clusters(i).PCA.weights = clusters(i).PCA.weights(1:size(clusters(i).all_data,1),:);
        %clusters(i).PCA.avgs = clusters(i).PCA.avgs(1:size(clusters(i).all_data,1),:);
    end

    % -------------------------------------------------------------------------------
    % for each group, separate into contiguous clusters
    % -------------------------------------------------------------------------------

    grps = unique(clusters(i).PCA.class(clusters(i).PCA.class~=0)); % values are component of origin
    for j = 1:length(grps),
        wh = find(clusters(i).PCA.class==grps(j));
        XYZ = clusters(i).XYZ(:,wh);
        cl_index = spm_clusters(XYZ) ./ 100;
        clusters(i).PCA.class(wh) = clusters(i).PCA.class(wh) + cl_index;
    end
    ngrps = length(unique(clusters(i).PCA.class));
    fprintf(1,'%3.0f contiguous clusters separated by class',ngrps)

    % -------------------------------------------------------------------------------
    % average within classes / contiguous regions
    % -------------------------------------------------------------------------------

    grps = unique(clusters(i).PCA.class(clusters(i).PCA.class~=0)); % values are component of origin
    for j = 1:length(grps),
        clusters(i).PCA.avgs(:,j) = mean(a(:,find(clusters(i).PCA.class==grps(j))),2);,
        freq(j) = sum(clusters(i).PCA.class == grps(j));
    end

    disp(['Cluster ' num2str(i) ', ' num2str(clusters(i).numVox) ' voxels: ' num2str(size(clusters(i).PCA.pcomps,2)) ' components'])
    fprintf(1,'\tMean\tEigval\tcorrel\t')

    % -------------------------------------------------------------------------------
    % * display each component and correlation with behavior
    % -------------------------------------------------------------------------------

    for j = 1:size(clusters(i).PCA.pcomps,2)

        %[H,P,CI,STATS] = TTEST(clusters(i).PCA.pcomps(:,j),0,.05,0);
        % skip the t-test.  t-tests on component scores don't make a lot of sense.
        fprintf(1,'\n\t%3.3f\t%3.3f\t',mean(clusters(i).PCA.pcomps(:,j)),clusters(i).PCA.eigval(j))

        if length(varargin) > 0
            if ~isempty(varargin{1})
                co = corrcoef(clusters(i).PCA.pcomps(:,j),varargin{1});
                co = co(1,2);
                fprintf(1,'%3.3f\t',co)
            end
        end

    end
    fprintf(1,'\n')

    % -------------------------------------------------------------------------------
    % * display classification info
    % -------------------------------------------------------------------------------
    disp(['Classified into ' num2str(ngrps) ' groups:'])
    fprintf(1,'\tClass\tVoxels\tMean\tt\tp\tcorrect. p\tcorrel\t')

    % for each component, test mean value and correlation with behavior
    for j = 1:length(grps)

        [H,P,CI,STATS] = TTEST(clusters(i).PCA.avgs(:,j),0,.05,0);
        fprintf(1,'\n\t%3.0f\t%3.0f\t%3.3f\t%3.3f\t%3.3f\t%3.3f\t',j,freq(j),mean(clusters(i).PCA.avgs(:,j)),STATS.tstat,P,P .* size(clusters(i).PCA.avgs,2))

        if length(varargin) > 0
            if ~isempty(varargin{1})
                co = corrcoef(clusters(i).PCA.avgs(:,j),varargin{1});
                co = co(1,2);
                fprintf(1,'%3.3f\t',co)
            end
        end

    end
    fprintf(1,'\n')

    % -------------------------------------------------------------------------------
    % * Plot, if requested
    % -------------------------------------------------------------------------------

    if plotflag,
        figure('Color','w'), subplot(1,3,1), imagesc(a), title(['Cl ' num2str(i) ': Data']), xlabel('Voxels'),ylabel('Subjects')
        subplot(1,3,2), imagesc(clusters(i).PCA.weights'), title(['Weights (eigenvectors)']), xlabel('Voxels'),ylabel('Eigenvectors')
        subplot(1,3,3), imagesc(clusters(i).PCA.pcomps), title(['Component scores (predictions)']), xlabel('Voxels'),ylabel('Subjects')

        a = [clusters(i).PCA.class' a']; a=sortrows(a,1); a = a(:,2:end)';
        figure;subplot 131; imagesc(a),title(['Cl ' num2str(i) ':Data sorted by class']), xlabel('Class'),ylabel('Subjects'),
        xlab = [sort(clusters(i).PCA.class(clusters(i).PCA.class~=0)) clusters(i).PCA.class(clusters(i).PCA.class==0)];
        set(gca,'XTick',1:length(clusters(i).PCA.class)); set(gca,'XTickLabel',xlab)
        subplot 132; imagesc(clusters(i).PCA.avgs),title('Class averages'), xlabel('Class'),ylabel('Subjects'),
        subplot 133; if length(varargin) > 0, if ~isempty(varargin{1}), imagesc(varargin{1}'), title('Behavior'),end,end

    end

    else
        disp(['Cluster ' num2str(i) ' has less than 3 voxels.'])
        clusters(i).PCA.class = ones(1,clusters(i).numVox);
        clusters(i).PCA.avgs = clusters(i).timeseries;
        grps = 1;
        if ~isfield(clusters,'correl'), clusters(1).correl = [];, end
    end

    % -------------------------------------------------------------------------------
    % * separate into subclusters, based on class membership
    % -------------------------------------------------------------------------------

    if ~isfield(clusters,'correl'), clusters(i).correl = [];, end
    disp('Recomputing correlations and z-scores and saving in subclusters')
    if ~isfield(clusters,'XYZ'), clusters(i).XYZ = ones(3,size(clusters(i).all_data,2));, end
    if ~isfield(clusters,'Z'), clusters(i).Z = ones(1,size(clusters(i).XYZ,2));, end
    clear subc

    for j = 1:length(grps)

        try
            subc(j) = clusters(i);
        catch
            warning('clusters does not have all required fields.'); break
        end

        wh = find(clusters(i).PCA.class == grps(j));

        if length(varargin) > 0
            if ~isempty(varargin{1})
                co = corrcoef(clusters(i).PCA.avgs(:,j),varargin{1});
                subc(j).correl = co(1,2);
            end
        end

        if isfield(subc(j),'title'), subc(j).title = [subc(j).title '_SUBCL_' num2str(j)];, end
        if isfield(subc(j),'name'), subc(j).name = [subc(j).name '_SUBCL_' num2str(j)];, end
        if isfield(subc(j),'numVox'), subc(j).numVox = length(wh);, end
        if isfield(subc(j),'Z'), subc(j).Z = subc(j).Z(wh);, end
        if isfield(subc(j),'XYZmm'), subc(j).XYZmm = subc(j).XYZmm(:,wh);, end
        if isfield(subc(j),'XYZ'), subc(j).XYZ = subc(j).XYZ(:,wh);, end

        if isfield(subc(j),'timeseries'), subc(j).timeseries = nanmean(subc(j).all_data(:,wh)')';, end
        %if isfield(subc(j),'snr'), subc(j).snr = subc(j).snr(wh);, end
        if isfield(subc(j),'center'), subc(j).center = center_of_mass(subc(j).XYZ,subc(j).Z);, end
        if isfield(subc(j),'mm_center'), subc(j).mm_center = center_of_mass(subc(j).XYZmm,subc(j).Z);, end

        if isfield(subc(j),'all_data'),
            subc(j).all_data = subc(j).all_data(:,wh);,
            for k = 1:size(subc(j).all_data,2)
                [H,P,CI,STATS] = ttest(subc(j).all_data(:,k),0,.05,0);
                subc(j).Z(k) = spm_t2z(STATS.tstat,STATS.df);
            end
        end

    end

    subclusters{i} = subc;

end

return

function [b,v,d,class] = pc(a,corrflag)
% a is original matrix, b is principal components, v is eigenvectors
% (weights on columns, which = weights on voxels)
% class is classification of voxels into groups based on component loadings

if corrflag, [v,d]=eig(corrcoef(a));, else, [v,d]=eig(cov(a));,end
b = (pinv(v) * a')' ./ repmat((diag(d)').^.5,size(a,1),1);
% i made this up: think of rptating each subject's scores (in cols of a')
% by the rotation matrix pinv(v), and normalizing by the sqrt of the eigenvalues
% pinv(v) and v are rotation matrices because det = 1, no shearing or dilation
%
% this appears to work to give scores as well
% both methods (above,below) are scaled versions of the splus factor scores
% the problem is that doing it two different ways in splus flips the signs
% of some components and not others (gui vs cmd line).

%X = a; R = corrcoef(a); A = v * (d^.5); B = inv(R) * A;
%scores = X * B;
% X is data, A is factor loading matrix, B is factor score coeff matrix
% this method, from the text, and the one giving b above produce identical results

A = v * (d^.5);

b = fliplr(b); v = fliplr(v); A = fliplr(A); %scores = fliplr(scores);

num = min(10,sum(diag(d) >= 1));
b = b(:,1:num); v = v(:,1:num); A = A(:,1:num);
origd = diag(d);
d = diag(d)'; d= fliplr(d); d = d(1:num);

if num == 0, warning('No eigenvalues above 1!');, origd, class = [];

else
    % classify each voxel into a group based on loading
    % use A, which re-introduces the comp variance, because we
    % want relationships with more variance to count more.
    % This just doesn't work so hot.  See docluster, below.

    %wh = A' == repmat(max(A'),size(A,2),1);
    %for i = 1:size(wh,2), tmp = find(wh(:,i));, class(i) = tmp(1); end
    %class(max(A') < .3) = 0;

end

%figure;plot(b,'r'),hold on;plot(a,'k'), hold on; plot(mean(a,2),'g--'),legend({'eig' 'orig' 'avg'})

return

function class = docluster(a,maxclusters,doplot)

    Y = pdist(a','euclid');     % transpose so the voxels are observations, eigenvectors the variables
    Z = linkage(Y,'complete');
    if maxclusters > 1
        class = cluster(Z,maxclusters)';

        if doplot,
            dendrogram(Z,0); title('Dendrogram for clustering')
        end

    else
        dendrogram(Z,0); title('Dendrogram for clustering')
        set(gcf,'Position',[10   601   800   500])
        maxclusters = input('Pick number of classes to save: ');

        if maxclusters == 1,
            class = ones(1,size(a,2));
        else
            class = cluster(Z,maxclusters)';
        end

    end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% SCRIPT TEST to view a manipulability ellipse

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
close all
fprintf('\nView manipulability ellipse and forces ellipse at q')

robot = load_robot('UR', 'UR10')

%q = [pi/2 pi/4 -pi/4 0.4 0.5 -pi/4]';
q = [0.15 0.15 0.15 0.15 0.15 0.15]';

%Plot manipulatiliby ellipse
drawrobot3d(robot, q)
hold on
%draw manipulability ellipse
plot_manipulability_ellipse(robot, q, 0.5)
%draw forces ellipse
%plot_forces_ellipse(robot, q, 0.1)

repo: ARTE


text:
function [innermost, inside] = find_innermost_boundary(bnd)

% FIND_INNERMOST_BOUNDARY locates innermost compartment of a BEM model
% by looking at the containment of the triangular meshes describing
% the surface boundaries
%
% [innermost] = find_innermost_boundary(bnd)
%
% with the boundaries described by a struct-array bnd with
%   bnd(i).pnt  vertices of boundary i (matrix of size Nx3)
%   bnd(i).tri  triangles of boundary i (matrix of size Mx3)

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

ncmp = length(bnd);

if ncmp==1
  innermost = 1;
  return
end

% try to locate the innermost compartment
for i=1:ncmp
for j=1:ncmp
  % determine for a single vertex on each surface if it is inside or outside the other surfaces
  curpos1 = bnd(i).pos(1,:); % any point on the boundary is ok
  curpos  = bnd(j).pos;
  curtri  = bnd(j).tri;
  if i==j
    inside(i,j) = 0;
  else
    inside(i,j) = bounding_mesh(curpos1, curpos, curtri);
  end
end
end
% assume that the sources are in the innermost compartment
tmp = sum(inside, 2);
[i, innermost] = max(tmp);


text:
% SHOW_EVENTS - Display events in epochs. Events selected by
%                 MAKE_TIMEWARP function can be optionally highlighted.
%                 Each epoch is visualized as a row in the output image with
%                 events marked by colored rectangles.
%
% Usage:
%       >> im = show_events(EEG, 'key1',value1,'key2',value2,...);
%
% Inputs:
%
%   EEG                 - dataset structure.
%
%
% Optional inputs (in 'key', value format):
%
%   'eventThicknessCoef'  - (positive number) adjust the thickness of event
%                           markers in the image. The default value is 1.
%                           Lower values reduce thickness and higher values
%                           increase it. For example, to decrease event
%                           rectangle width by half, set this parameter to
%                           0.5. {default: 1}
%   'eventNames'          - a cell array containing names of events to be
%                           displayed. Order and repetition are ignored.
%                           {default: all events if timeWarp is not provided)
%
%   'timeWarp'            - a structure with latencies (time-warp matrix
%                           for newtimef function) and epochs with correct
%                           sequence created by MAKE_TIMEWARP function.
%                           the subsets {default: false|0)
% Outputs:
%
%   im                    - color (RGB) image showing events of interest in
%                           epochs.
%
% Example:
%
%   % To display all events in data structure EEG
%   >> show_events(EEG);
%
%   % To highlight events selected by MAKE_TIMEWARP function with thin
%   %  event markers
%  >> show_events(EEG, 'eventThicknessCoef', 0.5, 'eventNames', timeWarp.eventSequence,
%   'timeWarp', timeWarp);
%
% Version 1.1
% Author: Nima Bigdely Shamlo, SCCN/INC/UCSD, 2008
% See also: MAKE_TIMEWARP, NEWTIMEF

% Copyright (C) Nima Bigdely Shamlo, SCCN/INC/UCSD, 2008
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function  im = show_events(EEG, varargin)

%% check inputs
EEG = change_events_to_string(EEG);

inputKeyValues = finputcheck(varargin, ...
   {'eventThicknessCoef'     'real'    [0 Inf]          1; ...
    'eventNames'     'cell'    {}        {}  ; ...
    'timeWarp'        'struct'    []          struct([]); ...
    });

eventThicknessCoef = 0;
eventNames = {};
timeWarp=[];

% place key values into function workspace variables
inputKeyValuesFields = fieldnames(inputKeyValues);
for i=1:length(inputKeyValuesFields)
    eval([inputKeyValuesFields{i} '= inputKeyValues.' inputKeyValuesFields{i} ';']);
end

if length(eventNames) == 0
    if length(timeWarp) == 0
        eventNames = [];
        for i=1:length(EEG.epoch)
            eventNames =  [eventNames EEG.epoch(i).eventtype];
        end
    else
        for i=1:length(timeWarp.eventSequence)
            if ischar(timeWarp.eventSequence{i})
                eventNames = [eventNames timeWarp.eventSequence{i}];
            else % in case it is a cell of strings
                for j=1:length(timeWarp.eventSequence{i})
                    eventNames = [eventNames timeWarp.eventSequence{i}{j}];
                end
            end
        end
        %        eventNames = timeWarp.eventSequence; % if event names is not provided only show timeWarp events
    end
end

if length(timeWarp) == 0
    timeWarp=[];
    timeWarp.latencies = [];
    timeWarp.epochs = 1:length(EEG.epoch);
end
%% set image parameters

imWidth = 300*4;
imHeight = 240*4;
im = double(zeros(imHeight, imWidth, 3));

%% find latencies for events of interest in each epoch

uniqueEventNames = unique_cell_string(eventNames);
for epochNumber = 1:length(EEG.epoch)
    for eventNumber = 1:length(uniqueEventNames)
       ids = find(strcmp(EEG.epoch(epochNumber).eventtype,uniqueEventNames{eventNumber})); % index for events with this type in the current epoch
       latency{epochNumber, eventNumber} = cell2mat(EEG.epoch(epochNumber).eventlatency(ids));
    end
end

%% find a good default value for event marker width based on quantiles
% of inter-event latencies.

intervals = [];
for i=1:size(latency,1)
    intervals = [intervals diff(sort(cell2mat(latency(i,:)),'ascend'))]; %#ok<AGROW>
end

if isempty(intervals)
    eventLineWidth = 0.05*eventThicknessCoef;
else
    eventLineWidth = round(eventThicknessCoef * 0.5 * imWidth * -quantile(-intervals,0.8) /  ((EEG.xmax - EEG.xmin)*1000));
end

%% place colored rectangle block o image to display events
eventColors = lines(length(uniqueEventNames));
for epochNumber = 1:length(EEG.epoch)
    for eventNumber = 1:length(uniqueEventNames)
        for eventInstanceNumber = 1:length(latency{epochNumber, eventNumber})
            startHeight =  max(1,round((epochNumber - 1) * imHeight/ length(EEG.epoch)));
            endHeight = min(imHeight, round((epochNumber) * imHeight/ length(EEG.epoch)));
            startWidth = max(1,round(-eventLineWidth/2 + imWidth * (latency{epochNumber, eventNumber}(eventInstanceNumber) - (EEG.xmin*1000))/ ((EEG.xmax - EEG.xmin)*1000)));
            endWidth =  min(imWidth,startWidth + eventLineWidth);
            colorForEvent =  eventColorBasedOnTypeAndAcceptance(eventColors, eventNumber, epochNumber, latency{epochNumber, eventNumber}(eventInstanceNumber), timeWarp);
            for c = 1:3 % to make sure the highlighted events are not overwritten by other events in the image, we use 'max' and give the brightest color priority
                im(startHeight:endHeight,startWidth:endWidth,c) = max(im(startHeight:endHeight,startWidth:endWidth,c), colorForEvent(c));
            end
        end
    end
end

%% plot image with legend and axis labels

figure;
hold on;
for i=1:length(uniqueEventNames)
    plot([0 0],rand * [0 0],'Color',eventColors(i,:),'linewidth',10);
    uniqueEventNames{i} = strrep(uniqueEventNames{i},'_','-'); % replace _ with - so it is displayed correctly in the legend
end

legend(uniqueEventNames, 'Location', 'NorthWest');

image(round([EEG.xmin*1000 EEG.xmax*1000 ]), [1 length(EEG.epoch)],im);
xlim(round([EEG.xmin*1000 EEG.xmax*1000 ]))
ylim([1 length(EEG.epoch)])
axis normal;
xlabel('Latency', 'fontsize',16);
xlabel('Latency', 'fontsize',16);
ylabel('Epochs', 'fontsize',16);

if nargout == 0  % suppress output if not requested
    clear im ;
end

function color = eventColorBasedOnTypeAndAcceptance(eventColors, eventNumber, epochNumber, eventLatency,timeWarp)
accepted = ismember_bc(epochNumber, timeWarp.epochs);
if ~isempty(timeWarp.latencies)
    matchedEpoch = find(timeWarp.epochs ==  epochNumber);
    accepted = accepted && ~isempty(find(timeWarp.latencies(matchedEpoch,:) == eventLatency, 1));
end

color = eventColors(eventNumber,:);

if ~accepted
    color = color*0.3; % dim the color of unaccepted events.
end

function EEG = change_events_to_string(EEG)
needChange = false;
for i=1:length(EEG.event)
    if ~ischar(EEG.event(i).type)
        EEG.event(i).type = num2str( EEG.event(i).type );
        needChange = true;
    end
end

if needChange
    for e=1:length(EEG.epoch)
        for i=1:length(EEG.epoch(e).eventtype)
            EEG.epoch(e).eventtype(i) = {num2str(cell2mat(EEG.epoch(e).eventtype(i)))};
        end
    end
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	Implemetation of the tracker described in paper
%	"MEEM: Robust Tracking via Multiple Experts using Entropy Minimization",
%   Jianming Zhang, Shugao Ma, Stan Sclaroff, ECCV, 2014
%
%	Copyright (C) 2014 Jianming Zhang
%
%	This program is free software: you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation, either version 3 of the License, or
%	(at your option) any later version.
%
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%
%	You should have received a copy of the GNU General Public License
%	along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%	If you have problems about this software, please contact: jmzhang@bu.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function resample(I_vf,step_size)
global sampler;
global svm_tracker;
global config;

if nargin < 2
    step_size = max(round(sampler.template_size(1:2)/5),1);
end

feature_map = imresize(I_vf,config.ratio,'nearest');
step_size = max(round(min(sampler.template_size(1:2))/4),1);
step_size = step_size([1 1]);
rect=svm_tracker.output;

upleft = round([rect(1)-sampler.roi(1)+1,rect(2)-sampler.roi(2)+1]);
if ~((upleft(1)<1) || (upleft(2)<1) || (round(upleft(1)+rect(3)-1)>size(I_vf,2)) || (round(upleft(2)+rect(4)-1)>size(I_vf,1)))
    sub_win=I_vf(round(upleft(2): (upleft(2)+rect(4)-1)),round(upleft(1): (upleft(1)+rect(3)-1)),:);
    output_feat = imresize(sub_win,config.template_sz);
    svm_tracker.output_feat = output_feat(:)';
else
    warning('tracking window outside of frame');
    keyboard
end

sampler.patterns_dt = [im2colstep(feature_map,sampler.template_size,[step_size, size(I_vf,3)])';...
    svm_tracker.output_feat];
temp = repmat(rect,[size(sampler.patterns_dt,1),1]);

[X Y] = meshgrid(1:step_size(2):size(feature_map,2)-sampler.template_size(2)+1,1:step_size(1):size(feature_map,1)-sampler.template_size(1)+1);
temp(1:end-1,1) = (X(:)-1)/config.ratio + sampler.roi(1);
temp(1:end-1,2) = (Y(:)-1)/config.ratio + sampler.roi(2);

%% compute cost table
left = max(round(temp(:,1)),round(rect(1)));
top = max(round(temp(:,2)),round(rect(2)));
right = min(round(temp(:,1)+temp(:,3)),round(rect(1)+rect(3)));
bottom = min(round(temp(:,2)+temp(:,4)),round(rect(2)+rect(4)));
ovlp = max(right - left,0).*max(bottom - top, 0);
sampler.costs = 1 - ovlp./(2*rect(3)*rect(4)-ovlp);

sampler.state_dt = temp;

end

author: flyers

repo: drone-tracking


text:
classdef Filter_Boundary < Filter
    methods (Static)
        function obj = create(settings)
            switch settings.filter
                case 'P1'
                    obj = Filter_P1_LevelSet();
                case 'PDE'
                    obj = Filter_PDE_LevelSet();
            end
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [I_hyper,I_rgb,I_gt,g,G,extra,wl,wl_ori,Sel,s,bbl,H1] = ...
    prepare_fusion_dataset(dataset, options)
%PREPARE_FUSION_DATASET Summary of this function goes here
%   Detailed explanation goes here
if nargin < 2
    options = [];
end

I_hyper = [];
I_rgb = [];

switch dataset
    case {'pavia_nonrigid_algo','pavia_nonrigid_gt','pavia_robustness'}
        if strcmp(dataset, 'pavia_nonrigid_gt')
            all_data = FusionDataLoader.PaviaFusion_GT_Reg();
        elseif strcmp(dataset, 'pavia_nonrigid_algo')
            all_data = FusionDataLoader.PaviaFusion_LSQfreeform_Reg();
        elseif strcmp(dataset, 'pavia_robustness')
            all_data = FusionDataLoader.PaviaFusion_robustness_LSQfreeform_Reg();
        end

        data_index = parse_param(options,'data_index',7);
        selected_data = all_data(data_index);

        I_gt = selected_data.I_gt;
        I_hyper = selected_data.I_hyper;
        wl = selected_data.wl;
        I_rgb = selected_data.I_rgb;
        s = selected_data.s; % s = (sx, sy)
        extra = selected_data.extra; % extra = [extra_x, extra_y]
        g = selected_data.g;
        G = selected_data.G;
        Sel = selected_data.S;
        H1 = selected_data.H1;
        bbl = true(1,length(wl));

        wl_ori = wl;
    case 'gulfport'
        load('gulfport_fusion.mat');
        I_rgb = double(I_rgb);
        wl_ori = wl;
    case 'salton_sea'
        load('salton_sea_roi.mat');
        I1 = imread('salton_sea_roi_3.2015.jpg');
        load('reg_result_salton_sea.mat');
        bbl = logical(bbl);

        wl_ori = wl;
        [I_hyper,wl,I_rgb,g,H1,extra,s,G,Sel] = interpret_reg_params(I,wl,bbl,...
            I1,U2,V2,rho,s2,T2,sigma2);

        I_gt = zeros(size(I_rgb,1), size(I_rgb,2), length(bbl));
    otherwise
end

H1 = H1/255;
I_rgb = I_rgb/255;

% the competing methods do not have the extra margin boundary
if extra(1) > 0
    I_rgb = remove_rgb_margin(I_rgb, extra);
    I_gt = remove_rgb_margin(I_gt, extra);
end

end

author: zhouyuanzxcv

repo: Hyperspectral