text:
% HOMOTRANS - homogeneous transformation of points
%
% Function to perform a transformation on homogeneous points/lines
% The resulting points are normalised to have a homogeneous scale of 1
%
% Usage:
%           t = homoTrans(P,v);
%
% Arguments:
%           P  - 3 x 3 or 4 x 4 transformation matrix
%           v  - 3 x n or 4 x n matrix of points/lines

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.

%  Peter Kovesi
%  School of Computer Science & Software Engineering
%  The University of Western Australia
%  pk @ csse uwa edu au
%  http://www.csse.uwa.edu.au/~pk
%
%  April 2000
%  September 2007

function t = homotrans(P,v);

    [dim,npts] = size(v);

    if ~all(size(P)==dim)
	error('Transformation matrix and point dimensions do not match');
    end

    t = P*v;  % Transform

    for r = 1:dim-1     %  Now normalise
	t(r,:) = t(r,:)./t(end,:);
    end

    t(end,:) = ones(1,npts);

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [ output_args ] = rcnn_assign_soft_cluster_score(expidx)

p = rcnn_exp_params(expidx);

fn = [p.expDir '/' p.shortName '/imdb/cache/clusters.mat'];
load(fn);

joint_no = 13;

vecs_j = vecs{joint_no};
vecs_j = vecs_j(:, 1:end-1);

C = clusters{joint_no};
k = size(C, 1);
num_all = size(vecs_j, 1);

%[idx, C] = spectralclustering(vecs_j, k);
%load('~/spectral_clustering', 'C', 'idx');

idx = zeros(num_all, 1);
for j = 1:num_all
    vec = vecs_j(j, :);
    diff = bsxfun(@minus, C, vec);
    [~, min_dist] = min(sum(diff.^2, 2));
    idx(j) = min_dist;
end

%assign soft scores
D = pdist(C);
min_dist = min(D);
sigma = min_dist/1.5;

figure(1);
hold on;

is_3d = size(vecs_j, 2) > 2;

% do PCA to display stuff
if is_3d
    [coeff,score,latent,tsquared,explained,mu] = pca(vecs_j);
    vecs_j = score;
end

for i = 1:k
    iset = (idx == i);
    if is_3d
        scatter3(vecs_j(iset, 1), vecs_j(iset, 2), vecs_j(iset, 3), 5, 'filled');
    else
        scatter(vecs_j(iset, 1), vecs_j(iset, 2), 5, 'filled');
    end
end

for i = 1:num_all
    idx = unidrnd(num_all);
    vec = vecs_j(idx, :);
    dists = pdist2(C,vec);
    scores = exp(-dists.^2/(sigma^2));
    scores = scores/sum(scores);
    if is_3d
        scatter3(vec(1), vec(2), vec(3), 20, [0 0 0], 'filled');
    else
        scatter(vec(1), vec(2), 20, [0 0 0], 'filled');
    end
    find(scores' > 0.1)
    pause;
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [out] = score_KS(mdv, hilo, lambda)
% Calculates KS score
%
% USAGE:
%
%    [out] = score_ridge(mdv, hilo, lambda, crossval)
%
% INPUTS:
%    mdv:         structure
%    hilo:        (0's and 1's), ideally there will be a similar # of each.
%
% OPTIONAL INPUTS:
%    lambda:      weighting, if the mean is less than lambda, the scores get weighted less, default = .02
%
% OUTPUT:
%    out:         score

if nargin < 3
    lambda = .02;
end

[nvars, npoints] = size(mdv);
loset = mdv(:,hilo==0);
hiset = mdv(:,hilo==1);

weights = std(mdv,0,2);

scores = zeros(nvars, 1);
parfor i = 1:nvars

    [h, p] = kstest2(loset(i,:), hiset(i,:));
    scores(i) = max(log(p), -708); % log(realmin) to avoid scores of inf.

end

scores2 = scores .*(1-exp(-weights/lambda));
out = -sum(scores2);

author: opencobra

repo: cobratoolbox


text:
classdef prtClassPlsda < prtClass
    % prtClassPlsda  Partial least squares discriminant classifier
    %
    %    CLASSIFIER = prtClassPlsda returns a Partial least squares
    %    discriminant classifier
    %
    %    CLASSIFIER = prtClassPlsda(PROPERTY1, VALUE1, ...) constructs a
    %    prtClassMAP object CLASSIFIER with properties as specified by
    %    PROPERTY/VALUE pairs.
    %
    %    A prtClassPlsda object inherits all properties from the abstract
    %    class prtClass. In addition is has the following properties:
    %
    %    nComponents  -  The number of components
    %    Bpls         -  The regression weights, estimated during training
    %    xMeans       -  The xMeans, estimated during training
    %    yMeans       -  The yMeana, estimated during training
    %
    %    trainingTechnique - Either 'simpls' or 'pls2' - the training
    %       technique to utilize.  See prtUtilSimpls and prtUtilPls2.
    %
    %    For information on the partial least squares discriminant
    %    algorithm, please refer to the following URL:
    %
    %    http://en.wikipedia.org/wiki/Partial_least_squares_regression
    %
    %    A prtClassPlsda object inherits the TRAIN, RUN, CROSSVALIDATE and
    %    KFOLDS methods from prtAction. It also inherits the PLOT method
    %    from prtClass.
    %
    %    Example:
    %
    %   TestDataSet = prtDataGenUnimodal;      % Create some test and
    %   TrainingDataSet = prtDataGenUnimodal;  % training data
    %   classifier = prtClassPlsda;           % Create a classifier
    %   classifier = classifier.train(TrainingDataSet);    % Train
    %   classified = run(classifier, TestDataSet);         % Test
    %   subplot(2,1,1);
    %   classifier.plot;
    %   subplot(2,1,2);
    %   [pf,pd] = prtScoreRoc(classified,TestDataSet);
    %   h = plot(pf,pd,'linewidth',3);
    %   title('ROC'); xlabel('Pf'); ylabel('Pd');
    %
    %    See also prtClass, prtClassLogisticDiscriminant, prtClassBagging,
    %    prtClassMap, prtClassCap, prtClassBinaryToMaryOneVsAll, prtClassDlrt,
    %    prtClassKnn, prtClassFld, prtClassRvm, prtClassGlrt,  prtClass

    % Copyright (c) 2013 New Folder Consulting
    %
    % Permission is hereby granted, free of charge, to any person obtaining a
    % copy of this software and associated documentation files (the
    % "Software"), to deal in the Software without restriction, including
    % without limitation the rights to use, copy, modify, merge, publish,
    % distribute, sublicense, and/or sell copies of the Software, and to permit
    % persons to whom the Software is furnished to do so, subject to the
    % following conditions:
    %
    % The above copyright notice and this permission notice shall be included
    % in all copies or substantial portions of the Software.
    %
    % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
    % OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
    % MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
    % NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
    % DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
    % OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
    % USE OR OTHER DEALINGS IN THE SOFTWARE.

    properties (SetAccess=private)
        name = 'Partial Least Squares Discriminant' % Partial Least Squares Discriminant
        nameAbbreviation = 'PLSDA' % PLSDA
        isNativeMary = true;  % True
    end

    properties
        % w is a DataSet.nDimensions x 1 vector of projection weights
        % learned during Fld.train(DataSet)
        nComponents = 2;
    end

    properties (SetAccess=protected)
        Bpls     % The prediction weights
        xScores % T
        yScores % U
        xVectors% P
        yVectors% Q

        yMeansFactor % Factor to be added into regression output (accounts for X means and yMeans);
        vipXY   % VIP scores including variations in X
        vipY    % VIP scores including only variations in Y
    end

    properties (Hidden)
        trainingTechnique = 'simpls'; %{'Simpls','pls2'};
        %         xMeans % Used for PLS2 (vs. SIMPLS)
    end

    methods

        function self = prtClassPlsda(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end

        function self = set.trainingTechnique(self,val)
            if ~any(strcmpi(val,{'pls2','simpls'}))
                error('prtClassPlsda:trainingTechnique','trainingTechnique must be one of {''pls2'',''simpls''}; string provided was: %s',val);
            end
            self.trainingTechnique = val;
        end
        function self = set.nComponents(self,val)
            if ~prtUtilIsPositiveInteger(val)
                error('prt:prtClassPlsda:nComponents','nComponents must be a positive integer');
            end
            self.nComponents = val;
        end

    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,DataSet)

            X = DataSet.getObservations;

            Y = DataSet.getTargetsAsBinaryMatrix;

            if DataSet.nClasses < 2
                warning('prt:prtClassPlda:unaryOrUnlabeled','Training dataset for PLSDA has %d classes. This may cause issues.',DataSet.nClasses)
            end

            maxComps = min(size(X));
            if self.nComponents > maxComps;
                self.nComponents = maxComps;
            end

            xMeans = mean(X,1);
            yMeans = mean(Y,1);
            X = bsxfun(@minus, X, xMeans);
            Y = bsxfun(@minus, Y, yMeans);
            switch self.trainingTechnique
                case 'simpls'
                    [self.Bpls, R, self.xVectors, self.yVectors, self.xScores, self.yScores] = prtUtilSimpls(X,Y,self.nComponents);

                case 'pls2'
                    [self.Bpls, W, P, Q, T, U, B] = prtUtilPls2(X,Y,self.nComponents);

                    self.yVectors = Q;
                    self.xVectors = P;
                    self.yScores = T*B;
                    self.xScores = T;

                otherwise
                    error('prtClassPlsda:trainingTechnique','Invalid trainingTechnique specified');
            end
            self.yMeansFactor = yMeans - xMeans*self.Bpls;

            ssT = diag(self.xScores'*self.yScores);
            ssT = ssT./sum(ssT(:));

            self.vipXY = sqrt(size(X,2)) * sum(bsxfun(@times, bsxfun(@rdivide,self.xVectors,sqrt(sum(self.xVectors.^2,2))).^2, ssT'),2);

            ssT = diag(self.yScores'*self.yScores);
            ssT = ssT./sum(ssT(:));

            self.vipY = sqrt(size(X,2)) * sum(bsxfun(@times, bsxfun(@rdivide,self.xVectors,sqrt(sum(self.xVectors.^2,2))).^2, ssT'),2);
        end

        function DataSet = runAction(self,DataSet)
            yOut = bsxfun(@plus,DataSet.getObservations*self.Bpls, self.yMeansFactor);
            DataSet = DataSet.setObservations(yOut);
        end

        function xOut = runActionFast(self,xIn,ds) %#ok<INUSD>
            xOut = bsxfun(@plus,xIn*self.Bpls, self.yMeansFactor);
        end
    end

    methods (Hidden)
        function str = exportSimpleText(self) %#ok<MANU>
            titleText = sprintf('%% prtClassPlsda\n');
            plsdBText = prtUtilMatrixToText(full(self.Bpls),'varName','plsdaWeights');
            plsdYText = prtUtilMatrixToText(full(self.yMeansFactor),'varName','yMeansFactor');
            str = sprintf('%s%s%s',titleText,plsdBText,plsdYText);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [ a, b ] = p11_lim ( )

%*****************************************************************************80
%
%% P11_LIM returns the integration limits for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%[2006]-"A GA-based feature selection and parameters optimization for
%support vector machines"

% (9/12/2020)

function GA = jGeneticAlgorithmTour(feat,label,opts)
% Parameters
CR        = 0.8;   % crossover rate
MR        = 0.01;  % mutation rate
Tour_size = 3;     % tournament size

if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'CR'), CR = opts.CR; end
if isfield(opts,'MR'), MR = opts.MR; end
if isfield(opts,'Ts'), Tour_size = opts.Ts; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = jInitialization(N,dim);
% Fitness
fit  = zeros(1,N);
fitG = inf;
for i = 1:N
  fit(i) = fun(feat,label,X(i,:),opts);
  % Best update
  if fit(i) < fitG
    fitG = fit(i);
    Xgb  = X(i,:);
  end
end
% Pre
curve = zeros(1,max_Iter);
curve(1) = fitG;
t = 2;
% Generations
while t <= max_Iter
  % Preparation
  Xc1   = zeros(1,dim);
  Xc2   = zeros(1,dim);
  fitC1 = ones(1,1);
  fitC2 = ones(1,1);
  z     = 1;
  for i = 1:N
    if rand() < CR
      % Select two parents
      k1 = jTournamentSelection(fit,Tour_size,N);
      k2 = jTournamentSelection(fit,Tour_size,N);
      % Store parents
      P1 = X(k1,:);
      P2 = X(k2,:);
      % Single point crossover
      ind = randi([1, dim - 1]);
      % Crossover between two parents
      Xc1(z,:) = [P1(1:ind), P2(ind + 1:dim)];
      Xc2(z,:) = [P2(1:ind), P1(ind + 1:dim)];
      % Mutation
      for d = 1:dim
        % First child
        if rand() < MR
          Xc1(z,d) = 1 - Xc1(z,d);
        end
        % Second child
        if rand() < MR
          Xc2(z,d) = 1 - Xc2(z,d);
        end
      end
      % Fitness
      fitC1(1,z) = fun(feat,label,Xc1(z,:),opts);
      fitC2(1,z) = fun(feat,label,Xc2(z,:),opts);
      z = z + 1;
    end
  end
  % Merge population
  XX = [X; Xc1; Xc2];
  FF = [fit, fitC1, fitC2];
  % Select N best solution
  [FF, idx] = sort(FF,'ascend');
  X         = XX(idx(1:N),:);
  fit       = FF(1:N);
  % Best agent
  if fit(1) < fitG
    fitG = fit(1);
    Xgb  = X(1,:);
  end
  % Save
  curve(t) = fitG;
  fprintf('\nGeneration %d Best (GA Tournament)= %f',t,curve(t))
  t = t + 1;
end
% Select features based on selected index
Pos   = 1:dim;
Sf    = Pos(Xgb == 1);
sFeat = feat(:,Sf);
% Store results
GA.sf = Sf;
GA.ff = sFeat;
GA.nf = length(Sf);
GA.c  = curve;
GA.f  = feat;
GA.l  = label;
end

%// Tournament Selection //
function Index = jTournamentSelection(fit,Tour_size,N)
% Random positions based on position & Tournament Size
Tour_idx  = randsample(N,Tour_size);
% Select ftiness value based on position selected by tournament
Tour_fit  = fit(Tour_idx);
% Get position of best ftiness value (win tournament)
[~, idx]  = min(Tour_fit);
% Store the position
Index     = Tour_idx(idx);
end

function X = jInitialization(N,dim)
% Initialize X vectors
X = zeros(N,dim);
for i = 1:N
  for d = 1:dim
    if rand() > 0.5
      X(i,d) = 1;
    end
  end
end
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function [rectx,recty,area,perimeter] = minboundrect(x,y,metric)
% minboundrect: Compute the minimal bounding rectangle of points in the plane
% usage: [rectx,recty,area,perimeter] = minboundrect(x,y,metric)
%
% arguments: (input)
%  x,y - vectors of points, describing points in the plane as
%        (x,y) pairs. x and y must be the same lengths.
%
%  metric - (OPTIONAL) - single letter character flag which
%        denotes the use of minimal area or perimeter as the
%        metric to be minimized. metric may be either 'a' or 'p',
%        capitalization is ignored. Any other contraction of 'area'
%        or 'perimeter' is also accepted.
%
%        DEFAULT: 'a'    ('area')
%
% arguments: (output)
%  rectx,recty - 5x1 vectors of points that define the minimal
%        bounding rectangle.
%
%  area - (scalar) area of the minimal rect itself.
%
%  perimeter - (scalar) perimeter of the minimal rect as found
%
%
% Note: For those individuals who would prefer the rect with minimum
% perimeter or area, careful testing convinces me that the minimum area
% rect was generally also the minimum perimeter rect on most problems
% (with one class of exceptions). This same testing appeared to verify my
% assumption that the minimum area rect must always contain at least
% one edge of the convex hull. The exception I refer to above is for
% problems when the convex hull is composed of only a few points,
% most likely exactly 3. Here one may see differences between the
% two metrics. My thanks to Roger Stafford for pointing out this
% class of counter-examples.
%
% Thanks are also due to Roger for pointing out a proof that the
% bounding rect must always contain an edge of the convex hull, in
% both the minimal perimeter and area cases.
%
%
% See also: minboundcircle, minboundtri, minboundsphere
%
%
% default for metric
if (nargin<3) || isempty(metric)
  metric = 'a';
elseif ~ischar(metric)
  error 'metric must be a character flag if it is supplied.'
else
  % check for 'a' or 'p'
  metric = lower(metric(:)');
  ind = strmatch(metric,{'area','perimeter'});
  if isempty(ind)
    error 'metric does not match either ''area'' or ''perimeter'''
  end

  % just keep the first letter.
  metric = metric(1);
end

% preprocess data
x=x(:);
y=y(:);

% not many error checks to worry about
n = length(x);
if n~=length(y)
  error 'x and y must be the same sizes'
end

% if var(x)==0

% start out with the convex hull of the points to
% reduce the problem dramatically. Note that any
% points in the interior of the convex hull are
% never needed, so we drop them.
if n>3

    %%%%%%%%%%%%%%%%%%%%%%%%%
    if (var(x)== 0|| var(y)==0)
        if var(x)== 0
            x = [x-1;x(1); x+1 ];
            y = [y ;y(1);y];
            flag = 1;
        else
            y = [y-1;y(1); y+1 ];
            x = [x ;x(1);x];
            flag = 1;
        end

    else
        flag = 0;
     %%%%%%%%%%%%%%%%%%%%%%
    edges = convhull(x,y);  % 'Pp' will silence the warnings

    end

  % exclude those points inside the hull as not relevant
  % also sorts the points into their convex hull as a
  % closed polygon

  %%%%%%%%%%%%%%%%%%%%
  if flag == 0
  %%%%%%%%%%%%%%%%%%%%

  x = x(edges);
  y = y(edges);
  %%%%%%%%%%%%%%%%%%
  end
  %%%%%%%%%%%%%
  % probably fewer points now, unless the points are fully convex
  nedges = length(x) - 1;
elseif n>1
  % n must be 2 or 3
  nedges = n;
  x(end+1) = x(1);
  y(end+1) = y(1);
else
  % n must be 0 or 1
  nedges = n;
end

% now we must find the bounding rectangle of those
% that remain.

% special case small numbers of points. If we trip any
% of these cases, then we are done, so return.
switch nedges
  case 0
    % empty begets empty
    rectx = [];
    recty = [];
    area = [];
    perimeter = [];
    return
  case 1
    % with one point, the rect is simple.
    rectx = repmat(x,1,5);
    recty = repmat(y,1,5);
    area = 0;
    perimeter = 0;
    return
  case 2
    % only two points. also simple.
    rectx = x([1 2 2 1 1]);
    recty = y([1 2 2 1 1]);
    area = 0;
    perimeter = 2*sqrt(diff(x).^2 + diff(y).^2);
    return
end
% 3 or more points.

% will need a 2x2 rotation matrix through an angle theta
Rmat = @(theta) [cos(theta) sin(theta);-sin(theta) cos(theta)];

% get the angle of each edge of the hull polygon.
ind = 1:(length(x)-1);
edgeangles = atan2(y(ind+1) - y(ind),x(ind+1) - x(ind));
% move the angle into the first quadrant.
edgeangles = unique(mod(edgeangles,pi/2));

% now just check each edge of the hull
nang = length(edgeangles);
area = inf;
perimeter = inf;
met = inf;
xy = [x,y];
for i = 1:nang
  % rotate the data through -theta
  rot = Rmat(-edgeangles(i));
  xyr = xy*rot;
  xymin = min(xyr,[],1);
  xymax = max(xyr,[],1);

  % The area is simple, as is the perimeter
  A_i = prod(xymax - xymin);
  P_i = 2*sum(xymax-xymin);

  if metric=='a'
    M_i = A_i;
  else
    M_i = P_i;
  end

  % new metric value for the current interval. Is it better?
  if M_i<met
    % keep this one
    met = M_i;
    area = A_i;
    perimeter = P_i;

    rect = [xymin;[xymax(1),xymin(2)];xymax;[xymin(1),xymax(2)];xymin];
    rect = rect*rot';
    rectx = rect(:,1);
    recty = rect(:,2);
  end
end
% get the final rect

% all done

end % mainline end

author: stupidZZ


text:
function varargout = pow2(varargin)

switch class(varargin{1})

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

    case 'char'

        operator = CreateBasicOperator('convex','increasing','positive','callback');
        operator.derivative = @(x)(log(2)*pow2(x));

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

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

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%   Distribution code Version 1.0 -- 02/31/2020 by Wei Liu Copyright 2020
%
%   The code is created based on the method described in the following paper
%   [1] "Real-time Image Smoothing via Iterative Least Squares", Wei Liu, Pingping Zhang,
%        Xiaolin Huang, Jie Yang, Chunhua Shen and Ian Reid, ACM Transactions on Graphics,
%        presented at SIGGRAPH 2020.
%
%   The code and the algorithm are for non-comercial use only.

%  ---------------------- Input------------------------
%  F:              input image, can be gray image or RGB color image
%  lambda:     \lambda in Eq.(1), control smoothing strength
%  p:              the power norm in the Charbonnier penalty in Eq. (2)
%  eps:          the small constant number in the Charbonnier penalty in Eq. (2)
%  iter:           iteration number of the ILS in Eq. (8)

%  ---------------------- Output------------------------
%  U:             smoothed image

function U =ILS_LNorm(F, lambda, p, eps, iter)

F = single(F); % 'single' precision is very important to reduce the computational cost

gamma = 0.5 * p - 1;
c =  p * eps^gamma;

[N, M, D] = size(F);
sizeI2D = [N, M];

otfFx = psf2otf_Dx(sizeI2D); % equal to otfFx = psf2otf(fx, sizeI2D) where fx = [1, -1];
otfFy = psf2otf_Dy(sizeI2D); % equal to otfFy = psf2otf(fy, sizeI2D) where fy = [1; -1];

Denormin = abs(otfFx).^2 + abs(otfFy ).^2;
Denormin = repmat(Denormin, [1, 1, D]);
Denormin = 1 + 0.5 * c * lambda * Denormin;

U = F; % smoothed image

Normin1 = fft2(U);

for k = 1: iter

    % Intermediate variables \mu update, in x-axis and y-axis direction
    u_h = [diff(U,1,2), U(:,1,:) - U(:,end,:)];
    u_v = [diff(U,1,1); U(1,:,:) - U(end,:,:)];

    mu_h = c .* u_h - p .* u_h .* (u_h .* u_h + eps) .^ gamma;
    mu_v = c .* u_v - p .* u_v .* (u_v .* u_v + eps) .^ gamma;

    % Update the smoothed image U
    Normin2_h = [mu_h(:,end,:) - mu_h(:, 1,:), - diff(mu_h,1,2)];
    Normin2_v = [mu_v(end,:,:) - mu_v(1, :,:); - diff(mu_v,1,1)];

    FU = (Normin1 + 0.5 * lambda * (fft2(Normin2_h + Normin2_v))) ./ Denormin;
    U = real(ifft2(FU));

    Normin1 = FU;  % This helps to further enlarge the smoothing strength

end

author: wliusjtu


text:
function [f,p,r] = compute_f(T,H)

  if length(T) ~= length(H),
    size(T)
    size(H)
  end;

  N = length(T);
  numT = 0;
  numH = 0;
  numI = 0;
  for n=1:N,
    Tn = (T(n+1:end))==T(n);
    Hn = (H(n+1:end))==H(n);
    numT = numT + sum(Tn);
    numH = numH + sum(Hn);
    numI = numI + sum(Tn .* Hn);
  end;
  p = 1;
  r = 1;
  f = 1;
  if numH > 0,
    p = numI / numH;
  end;
  if numT > 0,
    r = numI / numT;
  end;
  if (p+r) == 0,
    f = 0;
  else
    f = 2 * p * r / (p + r);
  end;

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [SPM,pKX] = spm_spm_local(SPM)
% [Re]ML Estimation of a General Linear Model
% FORMAT [SPM,pKX] = spm_spm_local(SPM)
%
% Michal and Brian spent a week comparing how SPM computes the beta values
% and comparing that with Rory's/mrVIsta's glm code.
%
% The structure sent in to spm_get_data() contains information that
% specifies how the data will be scaled upon return.  If you only send in
% the file name, the data are unscaled.  We did not determined, however,
% which entry of the structure causes spm_get_data() to scale the returned
% values.
%
% Given the same design matrix and data, the two calculations are
% identical. (We didn't check the covariance stuff, this is assuming no
% whitening).  This can be hard to see because of problems in reading the
% data with SPM.  The routine spm_get_data() apparently can scale the time
% series that it returns.  This has the consequence that the beta weights
% may differ by this scale factor.  But apart from this, the beta weights
% are the same in the two cases, when the design matrix is the same.
%
%         pKX is the pseudo-inverse of the and smoothed (by K) design
%         matrix.
%
% required fields of SPM:
%
% xY.VY - nScan x 1 struct array of mapped image volumes
%         Images must have the same orientation, voxel size and data type
%       - Any scaling should have already been applied via the image handle
%         scalefactors.
%
% xX    - Structure containing design matrix information
%       - Required fields are:
%         xX.X      - Design matrix (raw, not temporally smoothed)
%         xX.name   - cellstr of parameter names corresponding to columns
%                     of design matrix
%       - Optional fields are:
%         xX.K      - cell of session-specific structures (see spm_filter)
%                   - Design & data are pre-multiplied by K
%                     (K*Y = K*X*beta + K*e)
%                   - Note that K should not smooth across block boundaries
%                   - defaults to speye(size(xX.X,1))
%         xX.W      - Optional whitening/weighting matrix used to give
%                     weighted least squares estimates (WLS). If not specified
%                     spm_spm will set this to whiten the data and render
%                     the OLS estimtes maximum likelihood
%                     i.e. W*W' = inv(xVi.V).
%
% xVi   - structure describing intrinsic temporal non-sphericity
%       - required fields are:
%         xVi.Vi    - array of non-sphericity components
%                   - defaults to {speye(size(xX.X,1))} - i.i.d.
%                   - specifying a cell array of contraints (Qi)
%                     These contraints invoke spm_reml to estimate
%                     hyperparameters assuming V is constant over voxels.
%                     that provide a high precise estimate of xX.V
%       - Optional fields are:
%         xX.V      - Optional non-sphericity matrix.  Cov(e) = sigma^2*V
%                     If not specified spm_spm will compute this using
%                     a 1st pass to identify signifcant voxels over which
%                     to estimate V.  A 2nd pass is then used to re-estimate
%                     the parameters with WLS and save the ML estimates
%                     (unless xX.W is already specified)
%
% xM    - Structure containing masking information, or a simple column vector
%         of thresholds corresponding to the images in VY.
%       - If a structure, the required fields are:
%         xM.TH - nVar x nScan matrix of analysis thresholds, one per image
%         xM.I  - Implicit masking (0=>none, 1 => implicit zero/NaN mask)
%         xM.VM - struct array of mapped explicit mask image volumes
% 		- (empty if no explicit masks)
%               - Explicit mask images are >0 for valid voxels to assess.
%               - Mask images can have any orientation, voxel size or data
%                 type. They are interpolated using nearest neighbour
%                 interpolation to the voxel locations of the data Y.
%       - Note that voxels with constant data (i.e. the same value across
%         scans) are also automatically masked out.
%
%
% In addition, the global default UFp is used to set a critical
% F-threshold for selecting voxels over which the non-sphericity
% is estimated (if required)
%
%_______________________________________________________________________
%
% spm_spm is the heart of the SPM package. Given image files and a
% General Linear Model, it estimates the model parameters, variance
% hyerparameters, and smoothness of standardised residual fields, writing
% these out to disk in the current working directory for later
% interrogation in the results section. (NB: Existing analyses in the
% current working directory are overwritten).  This directory
% now becomes the working directory for this analysis and all saved
% images are relative to this directory.
%
% The model is expressed via the design matrix (xX.X). The basic model
% at each voxel is of the form is Y = X*B + e, for data Y, design
% matrix X, (unknown) parameters B, and residual errors e. The errors
% are assummed to have a normal distribution.
%
% Sometimes confounds (e.g. drift terms in fMRI) are necessary.  These
% can be specified directly in the design matrix or implicitly, in terms
% of a residual forming matrix K to give a generalised linear model
% K*Y = K*X*B + K*e.  In fact K can be any matrix (e.g. a convolution
% matrix).
%
% In some instances i.i.d. assumptions about errors do not hold.  For
% example, with serially correlated (fMRI) data or correlations among the
% levels of a factor in repeated measures designs.  This non-sphericity
% can be specified in terms of components (SPM.xVi.Vi{i}). If specified
% these covariance components will then be estimated with ReML (restricted
% maximum likelihood) hyperparameters.  This estimation assumes the same
% non-sphericity for voxels that exceed the global F-threshold. The ReML
% estimates can then used to whiten the data giving maximum likelihood (ML)
% or Gauss-Markov estimators.  This entails a second pass of the data
% with an augmented model K*W*Y = K*W*X*B + K*W*e where W*W' = inv(xVi.V).
% xVi.V is the non-sphericity based on the hyperparameter esimates.
% W is stored in xX.W and cov(K*W*e) in xX.V.  The covariance of the
% parameter estimates is then xX.Bcov = pinv(K*W*X)*xX.V*pinv(K*W*X)';
%
% If you do not want ML estimates but want to use ordinary least squares
% (OLS) then simply set SPM.xX.W to the identity matrix. Any non-sphericity
% V will still be estimated but will be used to adjust the degrees of freedom
% of the ensuing statistics using the Satterthwaite approximation (c.f.
% the Greenhouse-Giesser corrections)
%
% If [non-spherical] variance components Vi are not specified xVi.Vi and
% xVi.V default to the identity matrix (i.e. i.i.d). The parameters are
% then estimated by OLS.  In this instance the OLS and ML estimates are
% the same.
%
% Note that only a single voxel-specific hyperaprameter (i.e. variance
% component) is estimated, even if V is not i.i.d.  This means spm_spm
% always implements a fixed-effects model.
% Random effects models can be emulated using a multi-stage procedure:
% This entails summarising the data with contrasts such that the fixed
% effects in a second model on the summary data are those effects of
% interest (i.e. the population effects). This means contrasts are
% re-entered into spm_spm to make an inference (SPM) at the next
% level. At this higher hierarchial level the residual variance for the
% model contains the appropriate variance components from lower levels.
% See spm_RandFX.man for further details and below.
%
% Under the additional assumption that the standardised residual images
% are non-stationary standard Gaussian random fields, results from
% Random field theory can be applied to estimate the significance
% statistic images (SPM's) adjusting p values for the multiple tests
% at all voxels in the search volume. The parameters required for
% this random field correction are the volume, and Lambda, the covariance
% matrix of partial derivatives of the standardised error fields.
%
% spm_est_smoothness estimates the variances of the partial derivatives
% in the axis directions (the diagonal of Lambda). The covariances (off
% diagonal elements of Lambda) are assumed to be zero.
%
%                           ----------------
%
%
% The volume analsed is the intersection of the threshold masks,
% explicit masks and implicit masks. See spm_spm_ui for further details
% on masking options.
%
%
%-----------------------------------------------------------------------
%
% The output of spm_spm takes the form of an SPM.mat file of the analysis
% parameters, and 'float' flat-file images of the parameter and variance
% [hyperparameter] estimates. An 8bit zero-one mask image indicating the
% voxels assessed is also written out, with zero indicating voxels outside
% tha analysed volume.
%
%                           ----------------
%
% The following SPM.fields are set by spm_spm (unless specified)
%
%     xVi.V      - estimated non-sphericity trace(V) = rank(V)
%     xVi.h      - hyperparameters  xVi.V = xVi.h(1)*xVi.Vi{1} + ...
%     xVi.Cy     - spatially whitened <Y*Y'> (used by ReML to estimate h)
%     xVi.CY     - <(Y - <Y>)*(Y - <Y>)'>    (used by spm_spm_Bayes)
%
%                           ----------------
%
%     Vbeta     - struct array of beta image handles (relative)
%     VResMS    - file struct of ResMS image handle  (relative)
%     VM        - file struct of Mask  image handle  (relative)
%
%                           ----------------
%     xX.W      - if not specified W*W' = inv(x.Vi.V)
%     xX.V      - V matrix (K*W*Vi*W'*K') = correlations after K*W is applied
%     xX.xKXs   - space structure for K*W*X, the 'filtered and whitened'
%                 design matrix
%               - given as spm_sp('Set',xX.K*xX.W*xX.X) - see spm_sp
%     xX.pKX    - pseudoinverse of K*W*X, computed by spm_sp
%     xX.Bcov   - xX.pKX*xX.V*xX.pKX - variance-covariance matrix of
%                 parameter estimates
%		  (when multiplied by the voxel-specific hyperparameter ResMS)
%                 (of the parameter estimates. (ResSS/xX.trRV = ResMS)    )
%     xX.trRV   - trace of R*V, computed efficiently by spm_SpUtil
%     xX.trRVRV - trace of RVRV
%     xX.erdf   - effective residual degrees of freedom (trRV^2/trRVRV)
%     xX.nKX    - design matrix (xX.xKXs.X) scaled for display
%                 (see spm_DesMtx('sca',... for details)
%
%                           ----------------
%
%     xVol.M    - 4x4 voxel->mm transformation matrix
%     xVol.iM   - 4x4 mm->voxel transformation matrix
%     xVol.DIM  - image dimensions - column vector (in voxels)
%     xVol.XYZ  - 3 x S vector of in-mask voxel coordinates
%     xVol.S    - Lebesgue measure or volume       (in voxels)
%     xVol.R    - vector of resel counts           (in resels)
%     xVol.FWHM - Smoothness of components - FWHM, (in voxels)
%
%                           ----------------
%
% xCon  - See Christensen for details of F-contrasts.  These are specified
%         at the end of spm_spm after the non-sphericity V has been defined
%         or estimated. The fist contrast tests for all effects of interest
%         assuming xX.iB and xX.iG index confounding or nuisance effects.
%
%     xCon      - Contrast structure (created by spm_FcUtil.m)
%     xCon.name - Name of contrast
%     xCon.STAT - 'F', 'T' or 'P' - for F/T-contrast ('P' for PPMs)
%     xCon.c    - (F) Contrast weights
%     xCon.X0   - Reduced design matrix (spans design space under Ho)
%                 It is in the form of a matrix (spm99b) or the
%                 coordinates of this matrix in the orthogonal basis
%                 of xX.X defined in spm_sp.
%     xCon.iX0  - Indicates how contrast was specified:
%                 If by columns for reduced design matrix then iX0 contains the
%                 column indices. Otherwise, it's a string containing the
%                 spm_FcUtil 'Set' action: Usuall one of {'c','c+','X0'}
%                 (Usually this is the input argument F_iX0.)
%     xCon.X1o  - Remaining design space (orthogonal to X0).
%                 It is in the form of a matrix (spm99b) or the
%                 coordinates of this matrix in the orthogonal basis
%                 of xX.X defined in spm_sp.
%     xCon.eidf - Effective interest degrees of freedom (numerator df)
%     xCon.Vcon - ...for handle of contrast/ESS image (empty at this stage)
%     xCon.Vspm - ...for handle of SPM image (empty at this stage)
%
%                           ----------------
%
%
% The following images are written to file
%
% mask.{img,hdr}                                   - analysis mask image
% 8-bit (uint8) image of zero-s & one's indicating which voxels were
% included in the analysis. This mask image is the intersection of the
% explicit, implicit and threshold masks specified in the xM argument.
% The XYZ matrix contains the voxel coordinates of all voxels in the
% analysis mask. The mask image is included for reference, but is not
% explicitly used by the results section.
%
%                           ----------------
%
% beta_????.{img,hdr}                                 - parameter images
% These are 16-bit (float) images of the parameter estimates. The image
% files are numbered according to the corresponding column of the
% design matrix. Voxels outside the analysis mask (mask.img) are given
% value NaN.
%
%                           ----------------
%
% ResMS.{img,hdr}                    - estimated residual variance image
% This is a 32-bit (double) image of the residual variance estimate.
% Voxels outside the analysis mask are given value NaN.
%
%                           ----------------
%
% RVP.{img,hdr}                      - estimated resels per voxel image
% This is a 32-bit (double) image of the RESELs per voxel estimate.
% Voxels outside the analysis mask are given value 0.  These images
% reflect the nonstationary aspects the spatial autocorrelations.
%
%
%-----------------------------------------------------------------------
%
% References:
%
% Christensen R (1996) Plane Answers to Complex Questions
%       Springer Verlag
%
% Friston KJ, Holmes AP, Worsley KJ, Poline JP, Frith CD, Frackowiak RSJ (1995)
% ``Statistical Parametric Maps in Functional Imaging:
%   A General Linear Approach''
%       Human Brain Mapping 2:189-210
%
% Worsley KJ, Friston KJ (1995)
% ``Analysis of fMRI Time-Series Revisited - Again''
%       Neuroimage 2:173-181
%
%-----------------------------------------------------------------------
%
% For those interested, the analysis proceeds a "block" at a time,
% The block size conforms to maxMem that can be set as a global variable
% MAXMEM (in bytes) [default = 2^20]
%
%_______________________________________________________________________
% @(#)spm_spm.m	2.66 Andrew Holmes, Jean-Baptiste Poline, Karl Friston 03/03/27
SCCSid   = '2.66';

%-Say hello
%-----------------------------------------------------------------------
SPMid    = spm('FnBanner',mfilename,SCCSid);
Finter   = spm('FigName','Stats: estimation...'); spm('Pointer','Watch')

%-Get SPM.mat if necessary
%-----------------------------------------------------------------------
if nargin ==0
	swd     = spm_str_manip(spm_get(1,'SPM.mat','Select SPM.mat'),'H');
	load(fullfile(swd,'SPM.mat'));
	SPM.swd = swd;
end

%-Change to SPM.swd if specified
%-----------------------------------------------------------------------
try
	cd(SPM.swd);
end

%-Ensure data are assigned
%-----------------------------------------------------------------------
try
	SPM.xY.VY;
catch
	helpdlg({	'Please assign data to this design',...
			'Use fMRI under model specification'});
	spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
	return
end

%-Delete files from previous analyses
%-----------------------------------------------------------------------
if exist(fullfile('.','mask.img'),'file') == 2

	str   = {'Current directory contains SPM estimation files:',...
		 'pwd = ',pwd,...
		 'Existing results will be overwritten!'};

	abort = spm_input(str,1,'bd','stop|continue',[1,0],1);
	if abort
		spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
		return
	else
		str = sprintf('Overwriting old results\n\t (pwd = %s) ',pwd)
		warning(str)
		drawnow
	end
end

files = {	'mask.???','ResMS.???','RVP.???',...
		'beta_????.???','con_????.???','ResI_????.???',...
		'ess_????.???', 'spm?_????.???'};

for i=1:length(files)
	if any(files{i} == '*'|files{i} == '?' )
		[j,null] = spm_list_files(pwd,files{i});
		for i=1:size(j,1)
			spm_unlink(deblank(j(i,:)))
		end
	else
		spm_unlink(files{i})
	end
end

%=======================================================================
% - A N A L Y S I S   P R E L I M I N A R I E S
%=======================================================================

%-Initialise
%=======================================================================
fprintf('%-40s: %30s','Initialising parameters','...computing')    %-#
xX            = SPM.xX;
[nScan nBeta] = size(xX.X);

%-If xM is not a structure then assumme it's a vector of thresholds
%-----------------------------------------------------------------------
try
	xM = SPM.xM;
catch
	xM = -ones(nScan,1)/0;
end
if ~isstruct(xM)
	xM = struct(	'T',	[],...
			'TH',	xM,...
			'I',	0,...
			'VM',	{[]},...
			'xs',	struct('Masking','analysis threshold'));
end

%-Check confounds (xX.K) and non-sphericity (xVi)
%-----------------------------------------------------------------------
if ~isfield(xX,'K')
	xX.K  = 1;
end
try
	%-If covariance components are specified use them
	%---------------------------------------------------------------
	xVi   = SPM.xVi;
catch

	%-otherwise assume i.i.d.
	%---------------------------------------------------------------
	xVi   = struct(	'form',  'i.i.d.',...
			'V',	 speye(nScan,nScan));
end

%-Get non-sphericity V
%=======================================================================
try
	%-If xVi.V is specified proceed directly to parameter estimation
	%---------------------------------------------------------------
	V     = xVi.V;
	str   = 'parameter estimation';

catch
	% otherwise invoke ReML selecting voxels under i.i.d assumptions
	%---------------------------------------------------------------
	V     = speye(nScan,nScan);
	str   = '[hyper]parameter estimation';
end

%-Get whitening/Weighting matrix: If xX.W exists we will save WLS estimates
%-----------------------------------------------------------------------
try
	%-If W is specified, use it
	%-------------------------------------------------------
	W     = xX.W;   % imagesc(W); title('Whitening matrix')
catch
	if isfield(xVi,'V')

		% otherwise make W a whitening filter W*W' = inv(V)
		%-------------------------------------------------------
		[u s] = spm_svd(xVi.V);
		s     = spdiags(1./sqrt(diag(s)),0,nScan,nScan);
		W     = u*s*u';
		W     = W.*(abs(W) > 1e-6);
		xX.W  = sparse(W);
	else
		% unless xVi.V has not been estimated - requiring 2 passes
		%-------------------------------------------------------
		W     = speye(nScan,nScan);
		str   = 'hyperparameter estimation (1st pass)';
	end
end

%-Design space and projector matrix [pseudoinverse] for WLS
%=======================================================================

% The first 'Set' initializes the SVD of the design matrix properties. In
% our default analysis, W is the identity matrix. In our case, we set K so
% that it is an impulse.
xX.xKXs = spm_sp('Set',spm_filter(xX.K,W*xX.X));		% KWX

% Hence, it turns out the with our settings, KWX is just X
% figure(1); subplot(1,3,1), plot(spm_filter(xX.K,W*xX.X));
%            subplot(1,3,2), plot(W*xX.X)
%            subplot(1,3,3), plot(xX.xKXs.X)

% Compute the pseudo-inverse of xX.xKXs (from above)
xX.pKX  = spm_sp('x-',xX.xKXs);				% projector

% So, after all this nonsense, it is still the case that pKX remains
% the pseudoinverse of X.
% max(abs(xX.X(:) - xX.xKXs.X(:)))

%-If xVi.V is not defined compute Hsqr and F-threshold under i.i.d.
%-----------------------------------------------------------------------
if ~isfield(xVi,'V')
	Fcname = 'effects of interest';
	iX0    = [SPM.xX.iB SPM.xX.iG];
	xCon   = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs);
	X1o    = spm_FcUtil('X1o', xCon(1),xX.xKXs);
	Hsqr   = spm_FcUtil('Hsqr',xCon(1),xX.xKXs);
	trRV   = spm_SpUtil('trRV',xX.xKXs);
	trMV   = spm_SpUtil('trMV',X1o);
	UFp    = spm('GetGlobal','UFp');
	UF     = spm_invFcdf(1 - UFp,[trMV,trRV]);
end

%-Image dimensions and data
%=======================================================================
VY       = SPM.xY.VY;
M        = VY(1).mat;
DIM      = VY(1).dim(1:3)';
xdim     = DIM(1); ydim = DIM(2); zdim = DIM(3);
YNaNrep  = spm_type(VY(1).dim(4),'nanrep');

%-Maximum number of residual images for smoothness estimation
%-----------------------------------------------------------------------
global     defaults
MAXRES   = defaults.stats.maxres;

%-maxMem is the maximum amount of data processed at a time (bytes)
%-----------------------------------------------------------------------
MAXMEM   = defaults.stats.maxmem;
nSres    = min(nScan,MAXRES);
blksz    = ceil(MAXMEM/8/nScan);				%-block size
nbch     = ceil(xdim*ydim/blksz);				%-# blocks
nbch = 1;
blksz    = ceil(MAXMEM);				%-block size

fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...done')               %-#

%-Initialise output images (unless this is a 1st pass for ReML)
%=======================================================================
if isfield(xX,'W')
fprintf('%-40s: %30s','Output images','...initialising')     %-#

%-Intialise the name of the new mask : current mask & conditions on voxels
%-----------------------------------------------------------------------
VM    = struct(		'fname',	'mask.img',...
			'dim',		[DIM',spm_type('uint8')],...
			'mat',		M,...
			'pinfo',	[1 0 0]',...
			'descrip',	'spm_spm:resultant analysis mask');
VM    = spm_create_vol(VM);

%-Intialise beta image files
%-----------------------------------------------------------------------
Vbeta(1:nBeta) = deal(struct(...
			'fname',	[],...
			'dim',		[DIM',spm_type('float')],...
			'mat',		M,...
			'pinfo',	[1 0 0]',...
			'descrip',	''));
for i = 1:nBeta
	Vbeta(i).fname   = sprintf('beta_%04d.img',i);
	Vbeta(i).descrip = sprintf('spm_spm:beta (%04d) - %s',i,xX.name{i});
	spm_unlink(Vbeta(i).fname)
end
Vbeta = spm_create_vol(Vbeta,'noopen');

%-Intialise residual sum of squares image file
%-----------------------------------------------------------------------
VResMS = struct(	'fname',	'ResMS.img',...
			'dim',		[DIM',spm_type('double')],...
			'mat',		M,...
			'pinfo',	[1 0 0]',...
			'descrip',	'spm_spm:Residual sum-of-squares');
VResMS = spm_create_vol(VResMS,'noopen');

%-Intialise residual images
%-----------------------------------------------------------------------
VResI(1:nSres) = deal(struct(...
			'fname',	[],...
			'dim',		[DIM',spm_type('double')],...
			'mat',		M,...
			'pinfo',	[1 0 0]',...
			'descrip',	'spm_spm:Residual image'));

for i = 1:nSres
	VResI(i).fname   = sprintf('ResI_%04d.img', i);
	VResI(i).descrip = sprintf('spm_spm:ResI (%04d)', i);
	spm_unlink(VResI(i).fname);
end
VResI = spm_create_vol(VResI);
end % (xX,'W')

fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...initialised')        %-#

%=======================================================================
% - F I T   M O D E L   &   W R I T E   P A R A M E T E R    I M A G E S
%=======================================================================

%-Intialise variables used in the loop
%=======================================================================
xords = [1:xdim]'*ones(1,ydim); xords = xords(:)';  % plane X coordinates
yords = ones(xdim,1)*[1:ydim];  yords = yords(:)';  % plane Y coordinates
S     = 0;                                          % Volume (voxels)
s     = 0;                                          % Volume (voxels > UF)
Cy    = 0;					    % <Y*Y'> spatially whitened
CY    = 0;					    % <Y*Y'> for ReML
EY    = 0;					    % <Y>    for ReML
i_res = round(linspace(1,nScan,nSres))';	    % Indices for residual

%-Initialise XYZ matrix of in-mask voxel co-ordinates (real space)
%-----------------------------------------------------------------------
XYZ   = zeros(3,xdim*ydim*zdim);

%-Cycle over bunches blocks within planes to avoid memory problems
%=======================================================================
% spm_progress_bar('Init',100,str,'');

totalVox = 0;
for z = 1:zdim				%-loop over planes (2D or 3D data)

    % current plane-specific parameters
    %-------------------------------------------------------------------
    zords   = z*ones(xdim*ydim,1)';	%-plane Z coordinates
    CrBl    = [];			%-parameter estimates
    CrResI  = [];			%-normalized residuals
    CrResSS = [];			%-residual sum of squares
    Q       = [];			%-in mask indices for this plane

    for bch = 1:nbch			%-loop over blocks

	%-# Print progress information in command window
	%---------------------------------------------------------------
	str   = sprintf('Plane %3d/%-3d, block %3d/%-3d',z,zdim,bch,nbch);
	fprintf('\r%-40s: %30s',str,' ')                                %-#

	%-construct list of voxels in this block
	%---------------------------------------------------------------
	I     = [1:blksz] + (bch - 1)*blksz;		%-voxel indices
	I     = I(I <= xdim*ydim);			%-truncate
	xyz   = [xords(I); yords(I); zords(I)];		%-voxel coordinates
	nVox  = size(xyz,2);				%-number of voxels

	%-Get data & construct analysis mask
	%===============================================================
	fprintf('%s%30s',sprintf('\b')*ones(1,30),'...read & mask data')%-#
	Cm    = logical(ones(1,nVox));			%-current mask

	%-Compute explicit mask
	% (note that these may not have same orientations)
	%---------------------------------------------------------------
	for i = 1:length(xM.VM)

		%-Coordinates in mask image
		%-------------------------------------------------------
		j      = xM.VM(i).mat\M*[xyz;ones(1,nVox)];

		%-Load mask image within current mask & update mask
		%-------------------------------------------------------
		Cm(Cm) = spm_get_data(xM.VM(i),j(:,Cm)) > 0;
	end

	%-Get the data in mask, compute threshold & implicit masks
	%---------------------------------------------------------------
	Y     = zeros(nScan,nVox);
	for i = 1:nScan

		%-Load data in mask
		%-------------------------------------------------------
		if ~any(Cm), break, end			%-Break if empty mask
		Y(i,Cm)  = spm_get_data(VY(i),xyz(:,Cm));

        % tmp = Y(i,Cm); hist(tmp(:));

		Cm(Cm)   = Y(i,Cm) > xM.TH(i);		%-Threshold (& NaN) mask
		if xM.I & ~YNaNrep & xM.TH(i) < 0	%-Use implicit mask
			Cm(Cm) = abs(Y(i,Cm)) > eps;
        end
        max(Y(:))
	end

	%-Mask out voxels where data is constant
	%---------------------------------------------------------------
	Cm(Cm) = any(diff(Y(:,Cm),1));
	Y      = Y(:,Cm);				%-Data within mask
	CrS    = sum(Cm);				%-# current voxels

    figure(1); plot(Y);
    totalVox = CrS + totalVox
    title(sprintf('Current voxels TS %.0f, total: ',CrS,totalVox));
    pause(.2)
    foo = analyzeRead(VY(i).fname);
    hist(foo(:),100)
    imtool(foo/max(foo(:)))
	%===============================================================
	%-Proceed with General Linear Model (if there are voxels)
	%===============================================================
	if CrS

		%-Whiten/Weight data and remove filter confounds
		%-------------------------------------------------------
		fprintf('%s%30s',sprintf('\b')*ones(1,30),'filtering')	%-#

        % For us, W is the identity.
        % figure(1); imagesc(W)
        % In this case, KWY is just Y, the data
        %     max(abs(KWY(:) - Y(:)))
		KWY   = spm_filter(xX.K,W*Y);

		%-General linear model: Weighted least squares estimation
		%------------------------------------------------------
		fprintf('%s%30s',sprintf('\b')*ones(1,30),' estimation') %-#

        % Wandell notes:
        %
        %  SPM computes the beta weights using the formula
        %
        %    KWY= KX*beta
        %
        %  Here, pKX is the pseudo-inverse of KX. X is the design matrix,
        %  and K is the term that introduces the convolution with the HRF
        %
        %  KWY is the whitened, filtered, time series. I don't underwstand
        %  why K is on both sides.  Maybe, therefore, I don't understand K.
        %
        %  The formula our glm computes is Y = X*beta
        %  So the differences are that we included the HRF into X and they
        %  assume we did not.  The W is the identity, so it doesn't matter.
        %  If we set up the SPM so that K is the impulse, we should get the
        %  same.  If we set up the SPM so that we pass it X prior to
        %  filtering with the HRF, then we should also get the same.

        % So,
		beta  = xX.pKX*KWY;			%-Parameter estimates
        pKX   = xX.pKX;             % Wandell added for return

		res   = spm_sp('r',xX.xKXs,KWY);	%-Residuals
		ResSS = sum(res.^2);			%-Residual SSQ
		clear KWY				%-Clear to save memory

		%-If ReML hyperparameters are needed for xVi.V
		%-------------------------------------------------------
		if ~isfield(xVi,'V')

			%-F-threshold & accumulate spatially whitened Y*Y'
			%-----------------------------------------------
			j   = sum((Hsqr*beta).^2,1)/trMV > UF*ResSS/trRV;
			j   = find(j);
			if length(j)
				q  = size(j,2);
				s  = s + q;
				q  = spdiags(sqrt(trRV./ResSS(j)'),0,q,q);
				Y  = Y(:,j)*q;
				Cy = Cy + Y*Y';
			end

		end % (xVi,'V')

		%-if we are saving the WLS parameters
		%-------------------------------------------------------
		if isfield(xX,'W')

			%-sample covariance and mean of Y (all voxels)
			%-----------------------------------------------
			CY         = CY + Y*Y';
			EY         = EY + sum(Y,2);

			%-Save betas etc. for current plane as we go along
			%-----------------------------------------------
			CrBl 	   = [CrBl,    beta];
			CrResI     = [CrResI,  res(i_res,:)];
			CrResSS    = [CrResSS, ResSS];

		end % (xX,'W')
		clear Y				%-Clear to save memory

	end % (CrS)

	%-Append new inmask voxel locations and volumes
	%---------------------------------------------------------------
	XYZ(:,S + [1:CrS]) = xyz(:,Cm);		%-InMask XYZ voxel coords
	Q                  = [Q I(Cm)];		%-InMask XYZ voxel indices
	S                  = S + CrS;		%-Volume analysed (voxels)

    end   % (bch)

    %-Plane complete, write plane to image files (unless 1st pass)
    %===================================================================
    if isfield(xX,'W')

    	fprintf('%s%30s',sprintf('\b')*ones(1,30),'...saving plane')	%-#

	%-Write Mask image
	%-------------------------------------------------------------------
	j   = sparse(xdim,ydim);
	if length(Q), j(Q) = 1; end
	VM    = spm_write_plane(VM, j, z);

	%-Write beta images
	%-------------------------------------------------------------------
	j   = NaN*ones(xdim,ydim);
	for i = 1:nBeta
		if length(Q), j(Q) = CrBl(i,:); end
		Vbeta(i) = spm_write_plane(Vbeta(i), j, z);
	end

	%-Write residual images
	%-------------------------------------------------------------------
	for i = 1:nSres
		if length(Q), j(Q) = CrResI(i,:); end
		VResI(i) = spm_write_plane(VResI(i), j, z);
   	end

	%-Write ResSS into ResMS (variance) image scaled by tr(RV) above
	%-------------------------------------------------------------------
	if length(Q), j(Q) = CrResSS;    end
	VResMS  = spm_write_plane(VResMS,j,z);

    end % (xX,'W')

    %-Report progress
    %-------------------------------------------------------------------
    fprintf('%s%30s',sprintf('\b')*ones(1,30),'...done')             %-#
    spm_progress_bar('Set',100*(bch + nbch*(z - 1))/(nbch*zdim));

end % (for z = 1:zdim)
fprintf('\n')                                                        %-#
spm_progress_bar('Clear')

%=======================================================================
% - P O S T   E S T I M A T I O N   C L E A N U P
%=======================================================================
if S == 0, warning('No inmask voxels - empty analysis!'), end

%-close written image files (unless 1st pass)
%=======================================================================
if isfield(xX,'W')

	fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...closing files')  %-#
	VM              = spm_close_vol(VM);
	Vbeta           = spm_close_vol(Vbeta);
	VResI           = spm_close_vol(VResI);
	VResMS          = spm_close_vol(VResMS);

end % (xVi,'V')

%-average sample covariance and mean of Y (over voxels)
%-----------------------------------------------------------------------
CY          = CY/S;
EY          = EY/S;
CY          = CY - EY*EY';

%-If not defined, compute non-sphericity V using ReML Hyperparameters
%=======================================================================
if ~isfield(xVi,'V')

	%-REML estimate of residual correlations through hyperparameters (h)
	%---------------------------------------------------------------
	str    = 'Temporal non-sphericity (over voxels)';
	fprintf('%-40s: %30s\n',str,'...REML estimation') %-#
	Cy            = Cy/s;

	% ReML for separable designs and covariance components
	%---------------------------------------------------------------
	if isstruct(xX.K)
		m     = length(xVi.Vi);
		h     = zeros(m,1);
		V     = sparse(nScan,nScan);
		for i = 1:length(xX.K)

			% extract blocks from bases
			%-----------------------------------------------
			q     = xX.K(i).row;
			p     = [];
			Qp    = {};
			for j = 1:m
				if nnz(xVi.Vi{j}(q,q))
					Qp{end + 1} = xVi.Vi{j}(q,q);
					p           = [p j];
				end
			end

			% design space for ReML (with confounds in filter)
			%-----------------------------------------------
			Xp         = xX.X(q,:);
			try
				Xp = [Xp xX.K(i).X0];
			end

			% ReML
			%-----------------------------------------------
			fprintf('%-30s- %i\n','  ReML Block',i);
			[Vp,hp]  = spm_reml(Cy(q,q),Xp,Qp);
			V(q,q)   = V(q,q) + Vp;
			h(p)     = hp;
		end
	else
		[V,h] = spm_reml(Cy,xX.X,xVi.Vi);
	end

	% normalize non-sphericity and save hyperparameters
	%---------------------------------------------------------------
	V           = V*nScan/trace(V);
	xVi.h       = h;
	xVi.V       = V;			% Save non-sphericity xVi.V
	xVi.Cy      = Cy;			%-spatially whitened <Y*Y'>
	SPM.xVi     = xVi;			% non-sphericity structure

	% If xX.W is not specified use W*W' = inv(V) to give ML estimators
	%---------------------------------------------------------------
	if ~isfield(xX,'W')
		save SPM SPM
		clear
		load SPM
		SPM = spm_spm(SPM);
		return
	end
end

%-Use non-sphericity xVi.V to compute [effective] degrees of freedom
%=======================================================================
xX.V            = spm_filter(xX.K,spm_filter(xX.K,W*V*W')');	% KWVW'K'
[trRV trRVRV]   = spm_SpUtil('trRV',xX.xKXs,xX.V);		% trRV (for X)
xX.trRV         = trRV;						% <R'*y'*y*R>
xX.trRVRV       = trRVRV;					%-Satterthwaite
xX.erdf         = trRV^2/trRVRV;				% approximation
xX.Bcov         = xX.pKX*xX.V*xX.pKX';				% Cov(beta)

%-Set VResMS scalefactor as 1/trRV (raw voxel data is ResSS)
%-----------------------------------------------------------------------
VResMS.pinfo(1) = 1/xX.trRV;
VResMS          = spm_create_vol(VResMS,'noopen');

%-Create 1st contrast for 'effects of interest' (all if not specified)
%=======================================================================
Fcname          = 'effects of interest';
try
	iX0     = [xX.iB xX.iG];
catch
	iX0     = [];
end
xCon            = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs);

%-Append contrasts for fMRI - specified by SPM.Sess(s).Fc(i)
%-----------------------------------------------------------------------
if isfield(SPM,'Sess')
    for s = 1:length(SPM.Sess)
	for i = 1:length(SPM.Sess(s).Fc)
	    iX0           = 1:nBeta;
	    iX            = SPM.Sess(s).col(SPM.Sess(s).Fc(i).i);
	    iX0(iX)       = [];
	    Fcname        = sprintf('Sess(%d):%s',s,SPM.Sess(s).Fc(i).name);
	    xCon(end + 1) = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs);
	end
    end
end

%-Smoothness estimates of component fields and RESEL counts for volume
%=======================================================================
try
	FWHM = SPM.xVol.FWHM;
	VRpv = SPM.xVol.VRpv;
	R    = SPM.xVol.R;
catch
	[FWHM,VRpv] = spm_est_smoothness(VResI,VM);
	R           = spm_resels_vol(VM,FWHM)';
end

%-Delete the residuals images
%=======================================================================
for  i = 1:nSres,
	spm_unlink([spm_str_manip(VResI(i).fname,'r') '.img']);
	spm_unlink([spm_str_manip(VResI(i).fname,'r') '.hdr']);
	spm_unlink([spm_str_manip(VResI(i).fname,'r') '.mat']);
end

%-Compute scaled design matrix for display purposes
%-----------------------------------------------------------------------
xX.nKX        = spm_DesMtx('sca',xX.xKXs.X,xX.name);

fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...done')               %-#

%-Save remaining results files and analysis parameters
%=======================================================================
fprintf('%-40s: %30s','Saving results','...writing')                 %-#

%-place fields in SPM
%-----------------------------------------------------------------------
SPM.xVol.XYZ   = XYZ(:,1:S);			%-InMask XYZ coords (voxels)
SPM.xVol.M     = M;				%-voxels -> mm
SPM.xVol.iM    = inv(M);			%-mm -> voxels
SPM.xVol.DIM   = DIM;				%-image dimensions
SPM.xVol.FWHM  = FWHM;				%-Smoothness data
SPM.xVol.R     = R;				%-Resel counts
SPM.xVol.S     = S;				%-Volume (voxels)
SPM.xVol.VRpv  = VRpv;				%-Filehandle - Resels per voxel

SPM.Vbeta      = Vbeta;				%-Filehandle - Beta
SPM.VResMS     = VResMS;			%-Filehandle - Hyperparameter
SPM.VM         = VM;				%-Filehandle - Mask

SPM.xVi        = xVi;				% non-sphericity structure
SPM.xVi.CY     = CY;				%-<(Y - <Y>)*(Y - <Y>)'>

SPM.xX         = xX;				%-design structure
SPM.xM         = xM;				%-mask structure

SPM.xCon       = xCon;				%-contrast structure

SPM.SPMid      = SPMid;
SPM.swd        = pwd;

%-Save analysis parameters in SPM.mat file
%-----------------------------------------------------------------------
save SPM SPM

%=======================================================================
%- E N D: Cleanup GUI
%=======================================================================
fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...done')               %-#
spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
fprintf('%-40s: %30s\n','Completed',spm('time'))                     %-#
fprintf('...use the results section for assessment\n\n')             %-#

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [out_val,out_ind,nevals]=tt_stat(tt, varargin)
% Compute statistics (min,max) of a tt_tensor via maxvol cross
%   [out_val,out_ind,nevals]=tt_stat(tt, varargin)
%
% Tries to compute values and indices of the sought quantities, listed in
% varargin. Possible choices are (case insensitive):
%   'LM', 'SM': largest, resp. smallest in magnitude entries,
%   'LR', 'SR': entries having largest, resp. smallest real parts, or
%   'LI', 'SI': entries having largest, resp. smallest imaginary parts.
% out_val is a 1 x N vector of output values, listed in the same order as
% in varargin, out_ind is a d x N matrix of indices, where d is the
% dimension of tt. The number of inspected elements is returned in nevals.
%
% Additional parameters can be given in varargin in the form
% 'param_name', param_value:
%       'kickrank': enrich tt by random subspaces of sizes kickrank. Larger
%       kickrank may improve the accuracy in "difficult" cases, for the price
%       of larger neval. Default value is 0.
%       'nswp': number of alternating iterations (default 3)
%
% Examples:
%   Mimic the call to tt_abs_max, tt_max and tt_min simultaneously:
%       [out,ind,neval]=tt_stat(tt, 'lm', 'lr', 'sr');
%   The same with kickrank set to 5:
%       [out,ind,neval]=tt_stat(tt, 'lm', 'lr', 'sr', 'kickrank', 5);
%   Minimize the banana function on a grid 32768 x 32768:
%       n = 15;
%       x = tkron(tt_x(2,n)/2^n, tt_ones(2,n));
%       y = tkron(tt_ones(2,n), tt_x(2,n)/2^n);
%       f = (0.5-x).^2 + 100*(y-x.^2).^2;
%       v = tt_stat(f,'sm','nswp',11)
%   This should give v ~ 2e-13.
%
% The difference with tt_*max* is the usage of dnr^2 entries
% instead of dr^2 (outputs are sought over TT blocks, not density matrices).
% Besides, it can carry out several iterations with random enrichment.

kickrank = 0; % random enrich may actually perturb the values, but sometimes it converges faster
nswp = 3; % these seem to be enough quite often

if (numel(varargin)==0)
    error('at least one sought quantity is needed');
end;

% Distinguish soughts from parameters
soughts = cell(1,numel(varargin));
j = 1;
while (j<=numel(varargin))
    if    ((strcmpi(varargin{j}, 'lm'))||(strcmpi(varargin{j}, 'sm'))||...
           (strcmpi(varargin{j}, 'lr'))||(strcmpi(varargin{j}, 'sr'))||...
           (strcmpi(varargin{j}, 'li'))||(strcmpi(varargin{j}, 'si')))
       soughts{j} = varargin{j};
    end;
    if (strcmpi(varargin{j}, 'kickrank'))
        kickrank = varargin{j+1};
        j = j+1;
    end;
    if (strcmpi(varargin{j}, 'nswp'))
        nswp = varargin{j+1};
        j = j+1;
    end;
    j = j+1;
end;
soughts(cellfun('isempty', soughts))=[];

d = tt.d;
n = tt.n;
r = tt.r;
tt = core2cell(tt);
r0 = r; % r will become different below

ind_global = cell(d+1,1);
ind_local = cell(d+1,1); % for reducing units
out_val = nan*ones(1,numel(soughts));
out_ind = zeros(d,numel(soughts));
% reductions
XX0 = cell(d+1,1);
XX0{1} = 1; XX0{d+1} = 1;

% We need to subtract unitary vectors to exlude already considered indices
% and select fresh ones
% Prepare storages for them
units = cell(d,numel(soughts));
for j=1:numel(soughts)
    for i=1:d
        units{i,j} = zeros(n(i),1);
    end;
end;
XU = cell(d+1,numel(soughts));
for j=1:numel(soughts)
    XU{1,j}=1; XU{d+1,j}=1;
end;
nevals = 0;

% Start iteration
swp = 0; % Zeroth iteration will be just QR and MV
dir = -1;
i = d;
while (swp<=nswp)
    if (swp>0)
        % Compute statistics
        q = reshape(tt{i}, r0(i), n(i)*r0(i+1));
        q = XX0{i}*q;
        q = reshape(q, r(i)*n(i), r0(i+1));
        q = q*XX0{i+1};
        q = reshape(q, r(i)*n(i)*r(i+1), 1);
        nevals = nevals + r(i)*n(i)*r(i+1); % number of evaluations
        ind_new = zeros(numel(soughts), 3);
        % We need to add the found extrema indices to maxvol, since they
        % will likely give good estimates in the next dimensions.
        pos = 0;
        for j=1:numel(soughts)
            out_changed = false;
            if (strcmpi(soughts{j}, 'lm'))
                [val,ind]=max(abs(q));
                out_changed = (val>abs(out_val(j)))||(isnan(out_val(j)));
            end;
            if (strcmpi(soughts{j}, 'lr'))
                [val,ind]=max(real(q));
                out_changed = (val>real(out_val(j)))||(isnan(out_val(j)));
            end;
            if (strcmpi(soughts{j}, 'li'))
                [val,ind]=max(imag(q));
                out_changed = (val>imag(out_val(j)))||(isnan(out_val(j)));
            end;
            if (strcmpi(soughts{j}, 'sm'))
                [val,ind]=min(abs(q));
                out_changed = (val<abs(out_val(j)))||(isnan(out_val(j)));
            end;
            if (strcmpi(soughts{j}, 'sr'))
                [val,ind]=min(real(q));
                out_changed = (val<real(out_val(j)))||(isnan(out_val(j)));
            end;
            if (strcmpi(soughts{j}, 'si'))
                [val,ind]=min(imag(q));
                out_changed = (val<imag(out_val(j)))||(isnan(out_val(j)));
            end;
            if (out_changed)
                out_val(j) = q(ind);
                ind = tt_ind2sub([r(i), n(i), r(i+1)], ind);
                % Add the new index to the expansion set
                ind_new(pos+1,:) = ind;
                pos = pos+1;
                if (i>1)&&(i<d)
                    out_ind(:,j) = indexmerge(ind_global{i}(ind(1),:), ind(2), ind_global{i+1}(ind(3),:));
                elseif (i==1)
                    out_ind(:,j) = indexmerge(ind(2), ind_global{i+1}(ind(3),:));
                else % i==d
                    out_ind(:,j) = indexmerge(ind_global{i}(ind(1),:), ind(2));
                end;
                % Update units vectors
                for k=1:d
                    units{k,j}=zeros(n(k),1);
                    units{k,j}(out_ind(k,j))=1;
                end;
                % and their reductions
                for k=1:i-1
                    XU{k+1,j} = XU{k,j}*units{k,j}.';
                    XU{k+1,j} = reshape(XU{k+1,j}, r(k)*n(k), 1);
                    XU{k+1,j} = XU{k+1,j}(ind_local{k+1}, :);
                end;
                for k=d:-1:i+1
                    XU{k,j} = units{k,j}*XU{k+1,j};
                    XU{k,j} = reshape(XU{k,j}, 1, n(k)*r(k+1));
                    XU{k,j} = XU{k,j}(:, ind_local{k});
                end;
            end;
        end;
        ind_new = ind_new(1:pos, :);

        % Subtract units to get fresh indices
        for j=1:numel(soughts)
            qu = XU{i,j}*reshape(units{i,j}*XU{i+1,j}, 1, n(i)*r(i+1));
            qu = reshape(qu, r(i)*n(i)*r(i+1),1);
            q = q-qu*out_val(j);
        end;
    end;

    if (dir>0)&&(i<d)
        q = reshape(q, r(i)*n(i), r(i+1));
        rr = zeros(r(i)*n(i), pos);
        % Enrich it by the extrema indices (just add corresponding unit
        % vectors)
        if (pos>0)
            for j=1:pos
                rr(tt_sub2ind([r(i), n(i)], ind_new(j,1:2)), j) = 1;
            end;
        end;
        % Enrich it by random
        rr2 = randn(r(i)*n(i), kickrank);
        [q,~]=qr([q, rr, rr2],0);
        rnew = size(q,2);

        % Maxvol
        ind = maxvol2(q);
        qq = q(ind,:);
        q = q/qq;
        % Indices
        ind_local{i+1} = ind; % for units
        ind_global{i+1} = indexmerge(ind_global{i}, (1:n(i))');
        ind_global{i+1} = ind_global{i+1}(ind,:);
        r(i+1) = rnew;
        % Reduction
        XX0{i+1} = reshape(tt{i}, r0(i), n(i)*r0(i+1));
        XX0{i+1} = XX0{i}*XX0{i+1};
        XX0{i+1} = reshape(XX0{i+1}, r(i)*n(i), r0(i+1));
        XX0{i+1} = XX0{i+1}(ind,:);
        % For units
        for j=1:numel(soughts)
            XU{i+1,j} = XU{i,j}*reshape(units{i,j}, 1, n(i));
            XU{i+1,j} = reshape(XU{i+1,j}, r(i)*n(i), 1);
            XU{i+1,j} = XU{i+1,j}(ind, :);
        end;
    elseif (dir<0)&&(i>1)
        % QR and right indices
        if (swp==0)
            % Don't perturb the original solution, just QR
            q = reshape(tt{i}, r(i), n(i)*r(i+1));
            rr = zeros(n(i)*r(i+1), 0);
        else
            q = reshape(q, r(i), n(i)*r(i+1));
            % Enrich by the extrema indices
            rr = zeros(n(i)*r(i+1), pos);
            if (pos>0)
                for j=1:pos
                    rr(tt_sub2ind([n(i) r(i+1)], ind_new(j,2:3)), j) = 1;
                end;
            end;
        end;
        % Enrich by random
        rr2 = randn(n(i)*r(i+1), kickrank*double(swp>0));
        [q,rv]=qr([q.',rr,rr2], 0);
        rv = rv(:,1:r(i)).';
        q = q.';
        rnew = size(q,1);
        % Maxvol
        ind = maxvol2(q.');
        qq = q(:, ind);
        q = qq\q;
        rv = rv*qq;
        % Indices
        ind_local{i} = ind; % for units
        ind_global{i} = indexmerge((1:n(i))', ind_global{i+1});
        ind_global{i} = ind_global{i}(ind,:);
        if (swp==0)
            tt{i} = reshape(q, rnew, n(i), r(i+1));
            % Cast rv to the next block
            q = reshape(tt{i-1}, r(i-1)*n(i-1), r(i));
            tt{i-1} = reshape(q*rv, r(i-1), n(i-1), rnew);
            r0(i) = rnew;
        end;
        r(i) = rnew;
        % Reduction
        XX0{i} = reshape(tt{i}, r0(i)*n(i), r0(i+1));
        XX0{i} = XX0{i}*XX0{i+1};
        XX0{i} = reshape(XX0{i}, r0(i), n(i)*r(i+1));
        XX0{i} = XX0{i}(:,ind);
        for j=1:numel(soughts)
            XU{i,j} = units{i,j}*XU{i+1,j};
            XU{i,j} = reshape(XU{i,j}, 1, n(i)*r(i+1));
            XU{i,j} = XU{i,j}(:, ind);
        end;
    end;

    i = i+dir;
    if ((dir>0)&&(i==d))||((dir<0)&&(i==1))
        dir = -dir;
        swp = swp+1;
    end;
end;

end

function [J]=indexmerge(varargin)
% Merges two or three indices in the little-endian manner
sz1 = max(size(varargin{1},1),1);
sz2 = max(size(varargin{2},1),1);
sz3 = 1;
if (nargin>2) % Currently allows only 3
    sz3 = max(size(varargin{3}, 1), 1);
end;
% J1 goes to the fastest index, just copy it
J1 = repmat(varargin{1}, sz2*sz3, 1);
% J2 goes to the middle
J2 = reshape(varargin{2}, 1, []);
J2 = repmat(J2, sz1, 1); % now sz1 ones will be the fastest
J2 = reshape(J2, sz1*sz2, []);
J2 = repmat(J2, sz3, 1);
J = [J1,J2];
if (nargin>2)
    % J3 goes to the slowest
    J3 = reshape(varargin{3}, 1, []);
    J3 = repmat(J3, sz1*sz2, 1); % now sz1 ones will be the fastest
    J3 = reshape(J3, sz1*sz2*sz3, []);
    J = [J,J3];
end;
end

author: oseledets

repo: TT-Toolbox


text:
function C=tprod(A,B)
% Computes the tensor-product of two P dimensional tensors (A & B)
% All dimensions of the inputed tensor except for the 2nd must be the same.
% The algorithm follows from paper by ...
%
% INPUTS:
%
% A - n1 x n2 x n3 .... nP tensor
% B - n2 x n4 x n3 .... nP tensor
%
% OUTPUTS:
%
% C - n1 x n4 x n3 .... nP tensor
%
% Original author :  Misha Kilmer, Ning Hao
% Edited by       :  G. Ely, S. Aeron, Z. Zhang, ECE, Tufts Univ. 03/16/2015

sa = size(A);
sb = size(B);

faces = length(sa);

nfaces = prod(sb(3:end));

for i = 3:faces
    A = fft(A,[],i);
    B = fft(B,[],i);
end
sc = sa;
sc(2) = sb(2);
C = zeros(sc);
for i = 1:nfaces
    C(:,:,i) = A(:,:,i)*B(:,:,i);
end
for i = 3:faces
    C = ifft(C,[],i);
end

author: andrewssobral


text:
%  internal function

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)

function hdr = load_untouch_nii_hdr(fileprefix, machine, filetype)

   if filetype == 2
      fn = sprintf('%s.nii',fileprefix);

      if ~exist(fn)
         msg = sprintf('Cannot find file "%s.nii".', fileprefix);
         error(msg);
      end
   else
      fn = sprintf('%s.hdr',fileprefix);

      if ~exist(fn)
         msg = sprintf('Cannot find file "%s.hdr".', fileprefix);
         error(msg);
      end
   end

   fid = fopen(fn,'r',machine);

   if fid < 0,
      msg = sprintf('Cannot open file %s.',fn);
      error(msg);
   else
      fseek(fid,0,'bof');
      hdr = read_header(fid);
      fclose(fid);
   end

   return					% load_nii_hdr

%---------------------------------------------------------------------
function [ dsr ] = read_header(fid)

        %  Original header structures
	%  struct dsr
	%       {
	%       struct header_key hk;            /*   0 +  40       */
	%       struct image_dimension dime;     /*  40 + 108       */
	%       struct data_history hist;        /* 148 + 200       */
	%       };                               /* total= 348 bytes*/

    dsr.hk   = header_key(fid);
    dsr.dime = image_dimension(fid);
    dsr.hist = data_history(fid);

    %  For Analyze data format
    %
    if ~strcmp(dsr.hist.magic, 'n+1') && ~strcmp(dsr.hist.magic, 'ni1')
        dsr.hist.qform_code = 0;
        dsr.hist.sform_code = 0;
    end

    return					% read_header

%---------------------------------------------------------------------
function [ hk ] = header_key(fid)

    fseek(fid,0,'bof');

	%  Original header structures
	%  struct header_key                     /* header key      */
	%       {                                /* off + size      */
	%       int sizeof_hdr                   /*  0 +  4         */
	%       char data_type[10];              /*  4 + 10         */
	%       char db_name[18];                /* 14 + 18         */
	%       int extents;                     /* 32 +  4         */
	%       short int session_error;         /* 36 +  2         */
	%       char regular;                    /* 38 +  1         */
	%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
	%       };                               /* total=40 bytes  */
	%
	% int sizeof_header   Should be 348.
	% char regular        Must be 'r' to indicate that all images and
	%                     volumes are the same size.

    v6 = version;
    if str2num(v6(1))<6
       directchar = '*char';
    else
       directchar = 'uchar=>char';
    end

    hk.sizeof_hdr    = fread(fid, 1,'int32')';	% should be 348!
    hk.data_type     = deblank(fread(fid,10,directchar)');
    hk.db_name       = deblank(fread(fid,18,directchar)');
    hk.extents       = fread(fid, 1,'int32')';
    hk.session_error = fread(fid, 1,'int16')';
    hk.regular       = fread(fid, 1,directchar)';
    hk.dim_info      = fread(fid, 1,'uchar')';

    return					% header_key

%---------------------------------------------------------------------
function [ dime ] = image_dimension(fid)

	%  Original header structures
	%  struct image_dimension
	%       {                                /* off + size      */
	%       short int dim[8];                /* 0 + 16          */
        %       /*
        %           dim[0]      Number of dimensions in database; usually 4.
        %           dim[1]      Image X dimension;  number of *pixels* in an image row.
        %           dim[2]      Image Y dimension;  number of *pixel rows* in slice.
        %           dim[3]      Volume Z dimension; number of *slices* in a volume.
        %           dim[4]      Time points; number of volumes in database
        %       */
	%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
	%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
	%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
	%       short int intent_code;   % short int unused1;   /* 28 + 2 */
	%       short int datatype;              /* 30 + 2          */
	%       short int bitpix;                /* 32 + 2          */
	%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
	%       float pixdim[8];                 /* 36 + 32         */
	%	/*
	%		pixdim[] specifies the voxel dimensions:
	%		pixdim[1] - voxel width, mm
	%		pixdim[2] - voxel height, mm
	%		pixdim[3] - slice thickness, mm
	%		pixdim[4] - volume timing, in msec
	%					..etc
	%	*/
	%       float vox_offset;                /* 68 + 4          */
	%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
	%       float scl_inter;   % float funused1;      /* 76 + 4 */
	%       short slice_end;   % float funused2;      /* 80 + 2 */
	%       char slice_code;   % float funused2;      /* 82 + 1 */
	%       char xyzt_units;   % float funused2;      /* 83 + 1 */
	%       float cal_max;                   /* 84 + 4          */
	%       float cal_min;                   /* 88 + 4          */
	%       float slice_duration;   % int compressed; /* 92 + 4 */
	%       float toffset;   % int verified;          /* 96 + 4 */
	%       int glmax;                       /* 100 + 4         */
	%       int glmin;                       /* 104 + 4         */
	%       };                               /* total=108 bytes */

    dime.dim        = fread(fid,8,'int16')';
    dime.intent_p1  = fread(fid,1,'float32')';
    dime.intent_p2  = fread(fid,1,'float32')';
    dime.intent_p3  = fread(fid,1,'float32')';
    dime.intent_code = fread(fid,1,'int16')';
    dime.datatype   = fread(fid,1,'int16')';
    dime.bitpix     = fread(fid,1,'int16')';
    dime.slice_start = fread(fid,1,'int16')';
    dime.pixdim     = fread(fid,8,'float32')';
    dime.vox_offset = fread(fid,1,'float32')';
    dime.scl_slope  = fread(fid,1,'float32')';
    dime.scl_inter  = fread(fid,1,'float32')';
    dime.slice_end  = fread(fid,1,'int16')';
    dime.slice_code = fread(fid,1,'uchar')';
    dime.xyzt_units = fread(fid,1,'uchar')';
    dime.cal_max    = fread(fid,1,'float32')';
    dime.cal_min    = fread(fid,1,'float32')';
    dime.slice_duration = fread(fid,1,'float32')';
    dime.toffset    = fread(fid,1,'float32')';
    dime.glmax      = fread(fid,1,'int32')';
    dime.glmin      = fread(fid,1,'int32')';

    return					% image_dimension

%---------------------------------------------------------------------
function [ hist ] = data_history(fid)

	%  Original header structures
	%  struct data_history
	%       {                                /* off + size      */
	%       char descrip[80];                /* 0 + 80          */
	%       char aux_file[24];               /* 80 + 24         */
	%       short int qform_code;            /* 104 + 2         */
	%       short int sform_code;            /* 106 + 2         */
	%       float quatern_b;                 /* 108 + 4         */
	%       float quatern_c;                 /* 112 + 4         */
	%       float quatern_d;                 /* 116 + 4         */
	%       float qoffset_x;                 /* 120 + 4         */
	%       float qoffset_y;                 /* 124 + 4         */
	%       float qoffset_z;                 /* 128 + 4         */
	%       float srow_x[4];                 /* 132 + 16        */
	%       float srow_y[4];                 /* 148 + 16        */
	%       float srow_z[4];                 /* 164 + 16        */
	%       char intent_name[16];            /* 180 + 16        */
	%       char magic[4];   % int smin;     /* 196 + 4         */
	%       };                               /* total=200 bytes */

    v6 = version;
    if str2num(v6(1))<6
       directchar = '*char';
    else
       directchar = 'uchar=>char';
    end

    hist.descrip     = deblank(fread(fid,80,directchar)');
    hist.aux_file    = deblank(fread(fid,24,directchar)');
    hist.qform_code  = fread(fid,1,'int16')';
    hist.sform_code  = fread(fid,1,'int16')';
    hist.quatern_b   = fread(fid,1,'float32')';
    hist.quatern_c   = fread(fid,1,'float32')';
    hist.quatern_d   = fread(fid,1,'float32')';
    hist.qoffset_x   = fread(fid,1,'float32')';
    hist.qoffset_y   = fread(fid,1,'float32')';
    hist.qoffset_z   = fread(fid,1,'float32')';
    hist.srow_x      = fread(fid,4,'float32')';
    hist.srow_y      = fread(fid,4,'float32')';
    hist.srow_z      = fread(fid,4,'float32')';
    hist.intent_name = deblank(fread(fid,16,directchar)');
    hist.magic       = deblank(fread(fid,4,directchar)');

    return					% data_history

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
% The COBRAToolbox: testWritePajekNet.m
%
% Purpose:
%     - testwritePajekNet tests the writePajekNet function
%       that transforms a hypergraph into a graph (Pajek format)
%
% Author:
%     - Marouen BEN GUEBILA 09/02/2017

global CBTDIR

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testwritePajekNet'));
cd(fileDir);

model = getDistributedModel('ecoli_core_model.mat');

%test solver packages
solverPkgs = {'tomlab_cplex'};

for k = 1:length(solverPkgs)

    fprintf('   Testing writePajekNet using %s ... ', solverPkgs{k});

    % change the COBRA solver (LP)
    solverOK = changeCobraSolver(solverPkgs{k}, 'LP', 0);

    if solverOK == 1
        %load test data
        fileID = fopen('COBRAmodeltest.net', 'r');
        testData = fscanf(fileID, '%s');
        fclose(fileID);

        %call fucntion
        writePajekNet(model);

        %save produced data
        fileID = fopen('COBRAmodel.net', 'r');
        Data = fscanf(fileID, '%s');
        fclose(fileID);

        %compare with produced data
        assert(isequal(testData, Data));

        %delete data
        delete COBRAmodel.net
    end
end

% change the directory
cd(currentDir)


text:
clear all
cbmex('open');
cbmex('trialconfig',1,'absolute')

samplingFreq = 30000;
eventChannel = 1;
contChannel = 1;
EventLFPs = [];
EventLFPTotal = [];
contRange = -50:50;
pause(0.5);

while size(EventLFPTotal,1) < 500
    [spikeData, procTime, contData] = cbmex('trialdata',1);
    %Example channel to choose for perievent marks
    continuousData = contData{contChannel,3};
    spikeTimes = spikeData{eventChannel,2}- procTime * samplingFreq;

    if ~isempty(spikeTimes)
        eventSpikes = bsxfun(@plus, contRange, double(spikeTimes(1:end-3)));
        EventLFPs = continuousData(eventSpikes);
        EventLFPTotal = [EventLFPTotal; EventLFPs];
        EventLFPTotal = reshape(EventLFPTotal, size(eventSpikes));

    end
%    for i = 1:length(SpikeTimes)

%     	  try
%             EventLFPs(size(EventLFPs,1)+1,:) = ContinuousData(SpikeTimes(i):SpikeTimes(i)+100);
%         catch
%         end
%    end
    pause(0.5);
end
cbmex('close');

plot(mean(EventLFPTotal,1));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Mh = hmxHorzcat(varargin)
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : hmxHorzcat.m                                  |
%|    #    |   VERSION    : 0.42                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 15.07.2018                                    |
%| ( === ) |   SYNOPSIS   : H-Matrix horizontal concatenation             |
%|  `---'  |                                                              |
%+========================================================================+

%%% Input analysis
Ml = varargin{1};
Mr = varargin{2};
if (nargin == 2)
    if isequal(Ml.pos{1},Mr.pos{1})
        Ix = (1:size(Ml,1))';
        Iy = (1:size(Ml,2)+size(Mr,2))';
    else
        error('hmxHorzcat.m : unavailable case');
    end
else
    Ml = varargin{1};
    Mr = varargin{2};
    Ix = varargin{3};
    Iy = varargin{4};
end

%%% Preparation
% Particles
X = Ml.pos{1}(Ix,:);
Y = [Ml.pos{2} ; Mr.pos{2}];
Y = Y(Iy,:);

% Accuracy
tol = max(Ml.tol,Mr.tol);

% Initialiation
Mh = hmx(X,Y,tol);

% Admissibility
[isfar,Xdim,Ydim] = hmxFar(Mh);

% Y Indices
bool = (Iy<=size(Ml,2));
Iyl  = Iy(bool);
Iyr  = Iy(~bool)-size(Ml,2);

% Compression for far distances
if isfar
    % Compression
    [Al,Bl] = lowrank(Ml,Ix,Iyl);
    [Ar,Br] = lowrank(Mr,Ix,Iyr);

    % Concatenation
    A = [Al,Ar];
    B = zeros(size(A,2),length(Iy));
    B(1:size(Al,2),bool)      = Bl;
    B(size(Al,2)+1:end,~bool) = Br;

    % Recompression
    [A,B] = hmxQRSVD(A,B,Mh.tol);

    % Validation
    flag = (numel(A)+numel(B)) <= prod(size(Mh));

else
    flag = 0;
end

%%% Compression
if flag
    % Type
    Mh.typ = 1;

    % Low-rank
    Mh.dat = {A,B};

%%%% Full or sparse for smallest box (stopping criterion)
elseif (min(size(Mh)) < 100)
    % Type
    Mh.typ = 2;

    % Full
    Mh.dat = zeros(size(Mh));

    % Fill data
    Mh.dat(:,bool)  = full(Ml,Ix,Iyl);
    Mh.dat(:,~bool) = full(Mr,Ix,Iyr);

%%% H-Matrix (recursion)
else
    % Type
    Mh.typ = 0;

    % Subdivision for X
    [I1,I2] = hmxSubdivide(X,Xdim);
    Mh.row  = {I1 , I1 , I2 , I2};

    % Subdivision for Y
    [I1,I2] = hmxSubdivide(Y,Ydim);
    Mh.col  = {I1 , I2 , I1 , I2};

    % H-Matrix (recursion)
    for i = 1:4
        % Coordinates
        Ir = Mh.row{i};
        Ic = Mh.col{i};

        % Recursion
        Mh.chd{i} = hmxHorzcat(Ml,Mr,Ix(Ir),Iy(Ic));
    end

    % Fusion
    Mh = hmxFusion(Mh);
end
end

% function Mh = hmxHorzcat(Ml,Mr)
% %+========================================================================+
% %|                                                                        |
% %|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
% %|           openHmx is part of the GYPSILAB toolbox for Matlab           |
% %|                                                                        |
% %| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
% %| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
% %| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
% %| LICENCE   : This program is free software, distributed in the hope that|
% %| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
% %| redistribute and/or modify it under the terms of the GNU General Public|
% %| License, as published by the Free Software Foundation (version 3 or    |
% %| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
% %| is available, please contact us to activate a "pay for remove" option. |
% %| CONTACT   : matthieu.aussal@polytechnique.edu                          |
% %| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
% %|                                                                        |
% %| Please acknowledge the gypsilab toolbox in programs or publications in |
% %| which you use it.                                                      |
% %|________________________________________________________________________|
% %|   '&`   |                                                              |
% %|    #    |   FILE       : hmxHorzcat.m                                  |
% %|    #    |   VERSION    : 0.40                                          |
% %|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
% %|  ( # )  |   CREATION   : 14.03.2017                                    |
% %|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
% %| ( === ) |   SYNOPSIS   : Concatenate H-Matrix along 2nd dimension      |
% %|  `---'  |                with the rule                                 |
% %|         |                H-Matrix > Full > Compr                       |
% %+========================================================================+
%
% % Check dimensions
% if ~(isa(Ml,'hmx') && isa(Mr,'hmx'))
%     error('hmxHorzcat.m : unavailable case')
% elseif (size(Ml,1)~=size(Mr,1))
%     error('hmxHorzcat.m : Dimensions of matrices being concatenated are not consistent')
% elseif norm(Ml.pos{1}-Mr.pos{1},'inf') > 1e-12
%     error('hmxHorzcat.m : unavailable case')
% end
%
% % Initialization
% Mh = hmx(Ml.pos{1},[Ml.pos{2};Mr.pos{2}],Ml.tol);
%
%
% %%% [H-Matrix & H-Matrix] -> H-Matrix
% if (Ml.typ == 0) && (Mr.typ == 0)
%     % Particles distance average
%     Yl1 = mean(Ml.chd{1}.pos{2},1);
%     Yl2 = mean(Ml.chd{2}.pos{2},1);
%     Yr1 = mean(Mr.chd{1}.pos{2},1);
%     Yr2 = mean(Mr.chd{2}.pos{2},1);
%
%     % Concatenation pairs
%     D     = [norm(Yl1-Yr1) norm(Yl1-Yr2) norm(Yl2-Yr1) norm(Yl2-Yr2)];
%     [r,d] = min(D);
%     if (d == 1) || (d == 4)
%         I = [1 1 ; 2 2 ; 3 3 ; 4 4];
%     else
%         I = [1 2 ; 2 1 ; 3 4 ; 4 3];
%     end
%
%     % Construction (recursion)
%     Mh.typ = 0;
%     for i = 1:4
%         Mh.chd{i} = hmxHorzcat(Ml.chd{I(i,1)},Mr.chd{I(i,2)});
%         Mh.row{i} = Ml.row{i};
%         Mh.col{i} = [Ml.col{I(i,1)};size(Ml,2)+Mr.col{I(i,2)}];
%     end
%
% %     % Construction (recursion)
% %     Mh.typ = 0;
% %     for i = 1:4
% %         Mh.chd{i} = hmxHorzcat(Ml.chd{i},Mr.chd{i});
% %         Mh.row{i} = Ml.row{i};
% %         Mh.col{i} = [Ml.col{i};size(Ml,2)+Mr.col{i}];
% %     end
%
%
% %     % Split for too big full leaf
% %     for i = 1:4
% %         if (Mh.chd{i}.typ == 2) && (min(size(Mh.chd{i})) > 100)
% %             tmp       = Mh.chd{i}
% %             Mh.chd{i} = hmxBuilder(tmp.pos{1},tmp.pos{2},tmp.dat,tmp.tol);
% %         end
% %     end
% %
% %     % Fusion
% %     Mh = hmxFusion(Mh);
%
% % [H-Matrix & Compr] -> H-Matrix
% elseif (Ml.typ == 0) && (Mr.typ == 1)
%     Mr = hmxSplit(Ml.row,Mr);
%     Mh = hmxHorzcat(Ml,Mr);
%
% % [H-Matrix & Full] -> H-Matrix
% elseif (Ml.typ == 0) && (Mr.typ == 2)
%     Mr = hmxSplit(Ml.row,Mr);
%     Mh = hmxHorzcat(Ml,Mr);
%
%
% %%% [Compr & H-matrix] -> H-Matrix
% elseif (Ml.typ == 1) && (Mr.typ == 0)
%     Ml = hmxSplit(Mr.row,Ml);
%     Mh = hmxHorzcat(Ml,Mr);
%
% % [Compr & Compr] -> Compr
% elseif (Ml.typ == 1) && (Mr.typ == 1)
%     % Concatenation
%     nl = size(Ml.dat{1},2);
%     nr = size(Mr.dat{1},2);
%     A  = [Ml.dat{1} , Mr.dat{1}];
%     B  = [Ml.dat{2} , zeros(nl,size(Mr,2)) ;
%         zeros(nr,size(Ml,2)) , Mr.dat{2}] ;
%
%     % Recompression
%     [A,B] = hmxQRSVD(A,B,Ml.tol);
%
%     % Update
%     Mh.typ = 1;
%     Mh.dat = {A,B};
%
%
% % [Compr & Full] -> Full
% elseif (Ml.typ == 1) && (Mr.typ == 2)
%     Mh.typ = 2;
%     Mh.dat = [Ml.dat{1}*Ml.dat{2} , Mr.dat];
%
%
% %%% [Full & H-matrix] -> H-Matrix
% elseif (Ml.typ == 2) && (Mr.typ == 0)
%     Ml = hmxSplit(Mr.row,Ml);
%     Mh = hmxHorzcat(Ml,Mr);
%
% % [Full & Compr] -> Full
% elseif (Ml.typ == 2) && (Mr.typ == 1)
%     Mh.typ = 2;
%     Mh.dat = [Ml.dat , Mr.dat{1}*Mr.dat{2}];
%
% % [Full & Full] -> Full
% elseif (Ml.typ == 2) && (Mr.typ == 2)
%     Mh.typ = 2;
%     Mh.dat = [Ml.dat,Mr.dat];
%
% else
%     error('hmxHorzcat.m : unavailable case')
% end
% end
%
%
% function Mh = hmxSplit(row,Mh)
% % Subdivision for X
% Mh.row = row;
%
% % Subdivision for Y
% [~,~,Ydim] = hmxFar(Mh);
% [I1,I2]    = hmxSubdivide(Mh.pos{2},Ydim);
% Mh.col     = {I1 , I2 , I1 , I2};
%
% % H-Matrix (recursion)
% for i = 1:4
%     % Initialization
%     Mh.chd{i} = hmx(Mh.pos{1}(Mh.row{i},:),Mh.pos{2}(Mh.col{i},:),Mh.tol);
%
%     % H-Matrix
%     if (Mh.typ == 0)
%         error('hmxSplit.m : unavailable case')
%
%     % Compressed leaf
%     elseif (Mh.typ == 1)
%         % Subdivision
%         A = Mh.dat{1}(Mh.row{i},:);
%         B = Mh.dat{2}(:,Mh.col{i});
%
%         % Recompression
%         [A,B] = hmxQRSVD(A,B,Mh.tol);
%
%         % Update
%         Mh.chd{i}.typ = 1;
%         Mh.chd{i}.dat = {A,B};
%         Mh.chd{i}.tol = Mh.tol;
%
%     % Full leaf
%     elseif (Mh.typ == 2)
%         Mh.chd{i}.typ = 2;
%         Mh.chd{i}.dat = Mh.dat(Mh.row{i},Mh.col{i});
%         Mh.chd{i}.tol = Mh.tol;
%     end
% end
%
% % Erase data
% Mh.typ = 0;
% Mh.dat = [];
% end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [out] = snowfall_1(In,T,p1,varargin)
%snowfall_1

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

% Flux function
% ------------------
% Description:  Snowfall based on temperature threshold
% Constraints:  -
% @(Inputs):    p1   - temperature threshold below which snowfall occurs [oC]
%               T    - current temperature [oC]
%               In   - incoming precipitation flux [mm/d]
%               varargin(1) - smoothing variable r (default 0.01)

if size(varargin,2) == 0
    out = In.*(smoothThreshold_temperature_logistic(T,p1));
elseif size(varargin,2) == 1
    out = In.*(smoothThreshold_temperature_logistic(T,p1,varargin(1)));
end

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function eframes = vl_frame2oell(frames)
% VL_FRAMES2OELL   Convert a geometric frame to an oriented ellipse
%   EFRAME = VL_FRAME2OELL(FRAME) converts the generic FRAME to an
%   oriented ellipses EFRAME. FRAME and EFRAME can be matrices, with
%   one frame per column.
%
%   A frame is either a point, a disc, an oriented disc, an ellipse,
%   or an oriented ellipse. These are represented respectively by 2,
%   3, 4, 5 and 6 parameters each, as described in VL_PLOTFRAME().  An
%   oriented ellipse is the most general geometric frame; hence, there
%   is no loss of information in this conversion.
%
%   If FRAME is an oriented disc or ellipse, then the conversion is
%   immediate. If, however, FRAME is not oriented (it is either a
%   point or an unoriented disc or ellipse), then an orientation must
%   be assigned. The orientation is chosen in such a way that the
%   affine transformation that maps the standard oriented frame into
%   the output EFRAME does not rotate the Y axis. If frames represent
%   detected visual features, this convention corresponds to assume
%   that features are upright.
%
%   If FRAME is a point, then the output is an ellipse with null area.
%
%   See: feature frames,
%   VL_PLOTFRAME(), VL_HELP().

% Author: Andrea Vedaldi

% Copyright (C) 2013 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

[D,K] = size(frames) ;
eframes = zeros(6,K) ;

switch D
  case 2
    eframes(1:2,:) = frames(1:2,:) ;

  case 3
    eframes(1:2,:) = frames(1:2,:) ;
    eframes(3,:)   = frames(3,:) ;
    eframes(6,:)   = frames(3,:) ;

  case 4
    r = frames(3,:) ;
    c = r.*cos(frames(4,:)) ;
    s = r.*sin(frames(4,:)) ;

    eframes(1:2,:) = frames(1:2,:) ;
    eframes(3:6,:) = [c ; s ; -s ; c] ;

  case 5
    eframes(1:2,:) = frames(1:2,:) ;
    eframes(3:6,:) = mapFromS(frames(3:5,:)) ;

  case 6
    eframes = frames ;

  otherwise
     error('FRAMES format is unknown.') ;
end

% --------------------------------------------------------------------
function A = mapFromS(S)
% --------------------------------------------------------------------
% Returns the (stacking of the) 2x2 matrix A that maps the unit circle
% into the ellipses satisfying the equation x' inv(S) x = 1. Here S
% is a stacked covariance matrix, with elements S11, S12 and S22.
%
% The goal is to find A such that AA' = S. In order to let the Y
% direction unaffected (upright feature), the assumption is taht
% A = [a b ; 0 c]. Hence
%
%  AA' = [a^2, ab ; ab, b^2+c^2] = S.

A = zeros(4,size(S,2)) ;
a = sqrt(S(1,:));
b = S(2,:) ./ max(a, 1e-18) ;

A(1,:) = a ;
A(2,:) = b ;
A(4,:) = sqrt(max(S(3,:) - b.*b, 0)) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
% Copy this template configuration file to your VOT workspace.
% Enter the full path to the ECO repository root folder.

CFWCR_repo_path = 'your CFWCR path'

tracker_label = 'CFWCR';
tracker_command = generate_matlab_command('CFWCR_VOT()', {[CFWCR_repo_path, 'vot2017_trax']});
tracker_interpreter = 'matlab';
tracker_trax = true;

repo: CFWCR