text:
function H= plot_scalpEvolutionPlusChannel(erp, mnt, clab, ival, varargin)
%PLOT_SCALPEVOLUTIONPLUSCHANNEL - Display evolution of scalp topographies
%
%Synposis:
% H= plot_scalpEvolutionPlusChannel(ERP, MNT, CLAB, IVAL, <OPTS>)
%
%Description:
% Makes an ERP plot in the upper panel with given interval marked,
% and draws below scalp topographies for all marked intervals,
% separately for each Class. For each Classes topographies are plotted
% in one row and shared the same Color map scaling. (In future versions
% there might be different options for Color scaling.)
%
%Input:
% ERP: struct of epoched EEG data.
% MNT: struct defining an electrode montage
% CLAB: label of the channel(s) which are to be displayed in the
%       ERP plot.
% IVAL: [nIvals x 2]-sized array of interval, which are marked in the
%       ERP plot and for which scalp topographies are drawn.
%       When all interval are consequtive, ival can also be a
%       vector of interval borders.
% OPTS: struct or property/value list of optional fields/properties:
%  .IvalColor  - [nColors x 3]-sized array of rgb-coded Colors
%                with are used to mark intervals and corresponding
%                scalps. Colors are cycled, i.e., there need not be
%                as many Colors as interval. Two are enough,
%                default [0.4 1 1; 1 0.6 1].
%  .LegendPos  - specifies the position of the legend in the ERP plot,
%                default 0 (see help of legend for choices).
%  .PrintIvalUnits - appends the unit when writing the ival borders,
%                default 1
%
%The opts struct is passed to plot_scalpPattern
%
%Output:
% H: struct of handles to the created graphic objects.
%
%See also plot_scalpEvolution, plot_scalpPatternsPlusChannel, plot_scalp.

% Author(s): Benjamin Blankertz, Jan 2005

props= {'ChannelAtBottom', 0,                 'BOOL';
        'FigureColor',     [1 1 1],           'DOUBLE[3]';
        'GlobalCLim',      0,                 'BOOL';
        'IvalColor',       [.4 1 1; 1 .6 1],  'DOUBLE[- 3]';
        'LegendPos',       'Best',            'CHAR'
        'LineWidth',       3,                 'DOUBLE[1]';
        'PlotChannel',     1,                 'BOOL';
        'PrintIval',       0,                 'BOOL';
        'PrintIvalUnits',  1,                 'BOOL';
        'ScalePos',        'vert',            'CHAR';
        'ShrinkColorbar',  0,                 'DOUBLE';
        'Subplot',         [],                'DOUBLE|GRAPHICS';
        'SubplotChannel',  [],                'DOUBLE|GRAPHICS';
        'XUnit',           'ms',            'CHAR';
        'YUnit',           '\muV',          'CHAR';
};

props_scalpPattern= plot_scalpPattern;
props_channel= plot_channel;

if nargin==0,
  H= opt_catProps(props, props_scalpPattern, props_channel);
  return
end

% With input argument erp, we know it better
if util_getDataDimension(erp)==1
  props_channel= plotutil_channel1D;
else
  props_channel= plotutil_channel2D;
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props, props_scalpPattern, props_channel);

if isfield(erp, 'xUnit'),
  [opt,isdefault]= opt_overrideIfDefault(opt, isdefault, ...
                                         'XUnit', erp.xUnit);
end
if isfield(erp, 'yUnit'),
  [opt,isdefault]= opt_overrideIfDefault(opt, isdefault, ...
                                         'YUnit', erp.yUnit);
end

if isfield(opt, 'ColorOrder'),
  if isequal(opt.ColorOrder,'rainbow'),
    nChans= size(erp.y,1);
    opt.ColorOrder= hsv2rgb([(0.5:nChans)'/nChans ones(nChans,1)*[1 0.85]]);
%  else
%    if size(opt.ColorOrder,1)<size(erp.y,1),
%      opt.ColorOrder= repmat(opt.ColorOrder, [size(erp.y,1) 1]);
%    end
  end
else
  opt.ColorOrder= get(gca, 'ColorOrder');
end

opt_scalpPattern= opt_substruct(opt, props_scalpPattern(:,1));
opt_channel= opt_substruct(opt, props_channel(:,1));

[AxesStyle, dmy]= opt_extractPlotStyles(opt);

if max(sum(erp.y,2))>1,
  erp= proc_average(erp);
end

if size(ival,1)==1,
  ival= [ival(1:end-1)', ival(2:end)'];
end
if any(ival(:)>erp.t(end)),
  warning('interval out of epoch range: truncating');
  ival= min(ival, erp.t(end));
end
if any(ival(:)<erp.t(1)),
  warning('interval out of epoch range: truncating');
  ival= max(ival, erp.t(1));
end

nIvals= size(ival,1);
nColors= size(opt.IvalColor,1);
nClasses= length(erp.className);

mapCLim= zeros(2, nClasses, nIvals);
if opt.GlobalCLim,
  % determine common color range for all scalp maps
  commonCL= visutil_getCommonRange(erp, ival, 'CLim',opt_scalpPattern.CLim);
  mapCLim(1,:,:)= commonCL(1);
  mapCLim(2,:,:)= commonCL(2);
else
  if strcmp(opt.ScalePos, 'horiz'),
    % determine common color range for each interval
    for ii= 1:nIvals,
      commonCL= visutil_getCommonRange(erp, ival(ii,:), ...
                                       'CLim',opt_scalpPattern.CLim);
      mapCLim(1,:,ii)= commonCL(1);
      mapCLim(2,:,ii)= commonCL(2);
    end
  else
    % determine common color range for each class
    for cc= 1:nClasses,
      erp_cc= proc_selectClasses(erp, cc);
      commonCL= visutil_getCommonRange(erp_cc, ival, ...
                                       'CLim',opt_scalpPattern.CLim);
      mapCLim(1,cc,:)= commonCL(1);
      mapCLim(2,cc,:)= commonCL(2);
    end
  end
end

if isempty(opt.Subplot),
  clf;
end
set(gcf, 'Color',opt.FigureColor);

subplot_Offset= 0;
if opt.PlotChannel && ~isempty(clab),
  if ~isempty(opt.SubplotChannel),
    H.ax_erp= opt.SubplotChannel;
    axis_getQuietly(H.ax_erp);
  else
    if opt.ChannelAtBottom,
      H.ax_erp= subplotxl(1+nClasses, 1, 1+nClasses, ...
                          [0.1 0.05 0.1], [0.09 0 0.05]);
    else
      H.ax_erp= subplotxl(1+nClasses, 1, 1, 0.05, [0.09 0 0.05]);
      subplot_Offset= 1;
    end
  end
  if ~isempty(AxesStyle),
    set(H.ax_erp, AxesStyle{:});
  end
  hold on;   %% otherwise axis properties like ColorOrder are lost
  H.channel= plot_channel(erp, clab, opt_channel, 'Legend',0);
  for cc= 1:min(nColors,size(ival,1)),
    grid_markInterval(ival(cc:nColors:end,:), clab, opt.IvalColor(cc,:));
  end
  axis_redrawFrame(H.ax_erp);
  if ~isequal(opt.LegendPos, 'none'),
    if iscell(H.channel),
      hhh= H.channel{1};
    else
      hhh= H.channel;
    end
    % The following part can be removed. We do not allow numeric values
    % for LegendPos anymore.
    % Check matlab version for downward compatability
    if str2double(strtok(version,'.'))<7
      H.leg= legend(hhh.plot, erp.className, opt.LegendPos);
    else
      % check if LegendPos is integer
      if isnumeric(opt.LegendPos)
        switch(opt.LegendPos)
          case -1, loc = 'NorthEastOutside';
          case 0,  loc = 'Best';
          case 1,  loc = 'NorthEast';
          case 2,  loc = 'NorthWest';
          case 3,  loc = 'SouthWest';
          case 4,  loc = 'SouthEast';
          otherwise, loc = 'Best';
        end
        warning('Location numbers are obsolete, use "%s" instead of "%d"', ...
          loc,opt.LegendPos);
      else
        loc = opt.LegendPos;
      end
      H.leg= legend(hhh.plot, erp.className, 'Location',loc);
    end
  end
end

%if ~isempty(opt.Subplot),
%  opt.Subplot= reshape(opt.Subplot, [nClasses, nIvals]);
%end
cb_per_ival= strcmp(opt.ScalePos, 'horiz');
for cc= 1:nClasses,
  if ~any(any(erp.x(:,:,cc))),
    util_warning('empty_Class', sprintf('Class %d is empty', cc));
    continue;
  end
  for ii= 1:nIvals,
    if any(isnan(ival(ii,:))),
      continue;
    end
    if ~isempty(opt.Subplot),
      axis_getQuietly(opt.Subplot(cc, ii));
    else
      subplotxl(nClasses+opt.PlotChannel, nIvals, ...
                ii+(cc-1+subplot_Offset)*nIvals, ...
                [0.01+0.08*cb_per_ival 0.03 0.05], [0.05 0.02 0.1]);
    end
    opt_scalpPattern= setfield(opt_scalpPattern, 'ScalePos','none');
    opt_scalpPattern= setfield(opt_scalpPattern, 'Class',cc);
    H.scalp(cc,ii)= plot_scalpPattern(erp, mnt, ival(ii,:), ...
                           opt_scalpPattern, 'CLim',mapCLim(:,cc,ii));
    if cc==nClasses
      if opt.PrintIval,
        yLim= get(gca, 'yLim');
        if opt.PrintIvalUnits==2,
          ival_str= sprintf('%s - %s\n%s', num2str(ival(ii,1)), num2str(ival(ii,2)), opt.XUnit);
        elseif opt.PrintIvalUnits==1,
          ival_str= sprintf('%s - %s %s', num2str(ival(ii,1)), num2str(ival(ii,2)), opt.XUnit);
        else
          ival_str= sprintf('%s - %s', num2str(ival(ii,1)), num2str(ival(ii,2)));
        end;
        H.text_ival(ii)= text(mean(xlim), yLim(1)-0.04*diff(yLim), ival_str, ...
                              'verticalAli','top', 'horizontalAli','center');
      end
      if cb_per_ival,
        H.cb(ii)= plotutil_colorbarAside('horiz');
      end
    end
  end
  if strcmp(opt.ScalePos, 'vert'),
    H.cb(cc)= plotutil_colorbarAside;
    ylabel(H.cb(cc), ['[' opt.YUnit ']']);
    if opt.ShrinkColorbar>0,
      cbpos= get(H.cb(cc), 'Position');
      cbpos(2)= cbpos(2) + cbpos(4)*opt.ShrinkColorbar/2;
      cbpos(4)= cbpos(4) - cbpos(4)*opt.ShrinkColorbar;
      set(H.cb(cc), 'Position',cbpos);
    end
  end
  pos= get(H.scalp(cc,end).ax, 'position');
  yy= pos(2)+0.5*pos(4);
  H.background= visutil_getBackgroundAxis;
  H.text(cc)= text(0.01, yy, erp.className{cc});
  set(H.text(cc), 'verticalAli','top', ...
                  'horizontalAli','center', ...
                  'rotation',90, ...
                  'Visible','on', ...
                  'FontSize',12, ...
                  'fontWeight','bold');
  if isfield(opt, 'ColorOrder'),
    ccc= 1+mod(cc-1, size(opt.ColorOrder,1));
    set(H.text(cc), 'Color',opt.ColorOrder(ccc,:));
  end
end

if nargout<1,
  clear H
end

author: bbci


text:
 function x = ir_idct2(y)
%function x = ir_idct2(y)
%|
%| Compute 2D inverse DCT of many images
%| (matlab idct2 can handle only one image)
%|
%| in
%|	y [M N (many)]	2D (M by N) DCT coefficients for many images
%|
%| out
%|	x [M N (many)]	corresponding synthesized 2D images
%|
%| 2016 Anish Lahiri
%| 2016-08-024 modified by JF to include test and comments

if nargin < 1, ir_usage, end
if streq(y, 'test'), ir_idct2_test, return, end

[nx, ny, np] = size(y);
x = reshape(y, nx, []); % [M *]
x = idct(x); % [M *] 1D along 1st dim
x = permute(reshape(x, nx, ny, np), [2,1,3]); % [M N np]
x = reshape(x, ny, []); % [N *]
x = idct(x); % [N *] 1D along 2nd dim
x = permute(reshape(x, ny, nx, np), [2 1 3]); % [M N np]
x = reshape(x, size(y));

function ir_idct2_test
m = 8;
n = 10;
nrep = 100;
rng(0);
y = rand(m, n, nrep);
x1 = ir_idct2(y);
x2 = zeros(size(x1));
for ii=1:nrep
	x2(:,:,ii) = idct2(y(:,:,ii));
end
equivs(x1, x2)

y1 = ir_dct2(x1);
equivs(y, y1)

if im % show all the 2D basis functions
	n = 8;
	y = eye(n^2);
	y2 = reshape(y, n, n, []);
	x2 = ir_idct2(y2);
	im plc 1 2
	im(1, x2, '2D: correct')

	x1 = idct(y);
	x12 = reshape(x1, n, n, []);
	im(2, x12, '1D: incorrect')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
[~,structs]=libphaseretprotofile;

S = libstruct('phaseret_legla_init_params',struct());
calllib('libphaseret','phaseret_legla_init_params_defaults',S);

f = greasy;
a = 256;
M = 2*1024;
M2 = floor(M/2) + 1;
gl = M;
L = dgtlength(size(f,1),a,M);
g = firwin('blackman',gl);
gd = long2fir(gabdual(g,a,M),gl);
N = L/a;
maxit = 200;

corig = dgtreal(f,{'blackman',gl},a,M);
s = abs(corig) + 1i*zeros(size(corig));
cinitPtr = libpointer('doublePtr',complex2interleaved(s));

cout = zeros(2*M2,N);
coutPtr = libpointer('doublePtr',cout);

tic;
calllib('libphaseret','phaseret_legla_d',cinitPtr,g,L,gl,1,a,M,maxit,coutPtr);
toc;
 cout2 = interleaved2complex(coutPtr.Value);
 %cout2 = phaselockreal(cout2,a,M);

frec = idgtreal(cout2,{'dual',{'blackman',gl}},a,M);
s2 = dgtreal(frec,{'blackman',gl},a,M);
magnitudeerrdb(s,s2)

tic;
cout2 = legla(s,g,a,M,'modtrunc','onthefly','flegla','maxit',maxit);
toc;

frec = idgtreal(cout2,{'dual',{'blackman',gl}},a,M);

s2 = dgtreal(frec,{'blackman',gl},a,M);
magnitudeerrdb(s,s2)

clear S;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [b]=real(a)
%Compute a real part of TT-tensor
%   [B]=REAL(A)
%
% TT-Toolbox 2.2, 2009-2013
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%
% For details see Thm4 in Dolgov, Khoromskij, Savostyanov,
% "Superfast Fourier transform using the QTT format"
% http://dx.doi.org/10.1007/s00041-012-9227-4
%
%---------------------------

if (a.d==1)
    b=a;
    b.core=real(b.core);
    exit;
end
d=a.d;

%Determine the size and preallocate the result
b=tt_tensor;
b.r=a.r; b.r(2:d)=2*a.r(2:d);
b.d=a.d; b.n=a.n;
sz=dot(b.n.*b.r(1:d),b.r(2:d+1));
pos=(b.n.*b.r(1:d)).*b.r(2:d+1);
b.ps=cumsum([1;pos]);
b.core=zeros(sz,1);
apos=a.ps; bpos=b.ps;

% Fill cores
for i=1:d
    bb=zeros(b.r(i),b.n(i),b.r(i+1));
    aa=reshape(a.core(apos(i):apos(i+1)-1), [a.r(i),a.n(i),a.r(i+1)]);
    if (i==1)
        bb(1:a.r(i),:,1:a.r(i+1)           )=real(aa);
        bb(1:a.r(i),:,1+a.r(i+1):2*a.r(i+1))=imag(aa);
    elseif (i==d)
        bb(1:a.r(i),         :,1:a.r(i+1))= real(aa);
        bb(1+a.r(i):2*a.r(i),:,1:a.r(i+1))=-imag(aa);
    else
        bb(1:a.r(i),         :,1:a.r(i+1)           )= real(aa);
        bb(1:a.r(i),         :,1+a.r(i+1):2*a.r(i+1))= imag(aa);
        bb(1+a.r(i):2*a.r(i),:,1:a.r(i+1)           )=-imag(aa);
        bb(1+a.r(i):2*a.r(i),:,1+a.r(i+1):2*a.r(i+1))= real(aa);
    end
    b.core(bpos(i):bpos(i+1)-1)=bb(:);
end
return
end

author: oseledets

repo: TT-Toolbox


text:
function lambdas = hessianspectrum(problem, x, usepreconstr, storedb, key)
% Returns the eigenvalues of the (preconditioned) Hessian at x.
%
% function lambdas = hessianspectrum(problem, x)
% function lambdas = hessianspectrum(problem, x, useprecon)
% function lambdas = hessianspectrum(problem, x, useprecon, storedb)
% function lambdas = hessianspectrum(problem, x, useprecon, storedb, key)
%
% If useprecon is not set, or if it is set to 'noprecon' (default), this
% computes and returns the eigenvalues of the Hessian operator (which needs
% to be symmetric but not necessarily definite) on the tangent space at x.
% There are problem.M.dim() eigenvalues. Matlab's eigs is used internally.
%
% If useprecon is set to 'precon', the eigenvalues of the composition of
% the Hessian with the preconditioner at x are computed: Precon o Hessian.
% The preconditioner must have been defined in the problem structure and
% has to be symmetric, positive definite. It is supposed to approximate the
% inverse of the (Riemannian) Hessian. Ideally, the preconditioned Hessian
% is better conditioned (smaller ratio of largest to smallest eigenvalue in
% magnitude) than the non-preconditioned spectrum. The present tool can
% help assess that.
%
% The typical ways to define a preconditioner are via problem.precon or
% problem.sqrtprecon (see comment below). These should be function handles
% with the same input/output system as problem.hess for the Hessian.
%
% If the Hessian is not available from the problem structure, an
% approximate Hessian will be used. There are no guarantees of
% interpretability, but this may nevertheless be useful at times.
%
% Even though the Hessian and the preconditioner are both symmetric, their
% composition is not symmetric. This can slow down the call to 'eigs'
% substantially. If possible, you may specify the square root of the
% preconditioner in the problem structure, as sqrtprecon. This operator on
% the tangent space at x must also be symmetric, positive definite, and
% such that SqrtPrecon o SqrtPrecon = Precon. Then, the spectrum of the
% symmetric operator SqrtPrecon o Hessian o SqrtPrecon is computed: it is
% the same as the spectrum of Precon o Hessian, but is usually faster to
% compute. If both Precon and SqrtPrecon are provided, only SqrtPrecon will
% be used.
%
% The input and the output of the Hessian and of the preconditioner are
% projected on the tangent space to avoid undesired contributions of the
% ambient space.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% Requires the manifold description in problem.M to have these functions:
%
%   u_vec = vec(x, u_mat) :
%       Returns a column vector representation of the normal (usually
%       matrix) representation of the tangent vector u_mat. vec must be an
%       isometry between the tangent space (with its Riemannian metric) and
%       a subspace of R^n where n = length(u_vec), with the 2-norm on R^n.
%       In other words: it is an orthogonal projector.
%
%   u_mat = mat(x, u_vec) :
%       The inverse of vec (its adjoint).
%
%   u_mat_clean = tangent(x, u_mat) :
%       Subtracts from the tangent vector u_mat any component that would
%       make it "not really tangent", by projection.
%
%   answer = vecmatareisometries() :
%       Returns true if the linear maps encoded by vec and mat are
%       isometries, false otherwise. It is better if the answer is yes.
%
% See also: hessianextreme canGetPrecon canGetSqrtPrecon

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, July 3, 2013.
% Contributors:
% Change log:
%
%   Dec. 18, 2014 (NB):
%       The lambdas are now sorted when they are returned.
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%       Does no longer accept sqrtprecon as an input: the square root of
%       the preconditioner may now be specified directly in the problem
%       structure, following the same syntax as the preconditioner precon.
%
%   April 4, 2015 (NB):
%       By default, the spectrum is computed without the preconditioner's
%       effect, even if it is available. A new input option allows to
%       switch this behavior without the need to change the problem
%       structure.

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    % Manage the option to use or not use a preconditioner.
    % The input is a string. It is here transformed into a Boolean.
    if ~exist('usepreconstr', 'var') || isempty(usepreconstr)
        usepreconstr = 'noprecon';
    end
    switch lower(usepreconstr)
        case 'noprecon'
            useprecon = false;
        case 'precon'
            useprecon = true;
        otherwise
            % A bit of legacy code heads up.
            if isa(usepreconstr, 'function_handle')
                warning('manopt:hessianspectrum:oldsyntax', ...
                        ['This function no longer expects sqrtprecon ' ...
                         'as input. Place it in the problem structure.']);
            end
            error('Input useprecon must be either ''precon'' or ''noprecon''.');
    end

    % No warning if an approximate Hessian is available, as then the user
    % is presumably aware of what they are doing.
    if ~canGetHessian(problem) && ~canGetApproxHessian(problem)
        warning('manopt:hessianspectrum:nohessian', ...
                ['The Hessian appears to be unavailable.\n' ...
                 'Will try to use an approximate Hessian instead.\n'...
                 'Since this approximation may not be linear or '...
                 'symmetric,\nthe computation might fail and the '...
                 'results (if any)\nmight make no sense.']);
    end

    vec = @(u_mat) problem.M.vec(x, u_mat);
    mat = @(u_vec) problem.M.mat(x, u_vec);
    tgt = @(u_mat) problem.M.tangent(x, u_mat);

    % n: size of a vectorized tangent vector
    % dim: dimension of the tangent space
    % necessarily, n >= dim.
    % The vectorized operators we build below will have at least n - dim
    % zero eigenvalues.
    n = length(vec(problem.M.zerovec(x)));
    dim = problem.M.dim();

    % It is usually a good idea to force a gradient computation to make
    % sure precomputable things are precomputed.
    if canGetGradient(problem)
        [unused1, unused2] = getCostGrad(problem, x, storedb, key); %#ok
    end

    hess = @(u_mat) tgt(getHessian(problem, x, tgt(u_mat), storedb, key));
    hess_vec = @(u_vec) vec(hess(mat(u_vec)));

    % Regardless of preconditioning, we can only have a symmetric
    % eigenvalue problem if the vec/mat pair of the manifold is an
    % isometry:
    vec_mat_are_isometries = problem.M.vecmatareisometries();

    if ~useprecon

        % No preconditioner to use: simply use the Hessian as is.

        eigs_opts.issym = vec_mat_are_isometries;
        eigs_opts.isreal = true;
        lambdas = eigs(hess_vec, n, dim, 'LM', eigs_opts);

    elseif canGetSqrtPrecon(problem)

        % There is a preconditioner, and we have its square root: deal with
        % the symmetric composition SqrtPrecon o Hessian o SqrtPrecon.

        sqrtprec = @(u_mat) tgt(getSqrtPrecon(problem, x, tgt(u_mat), storedb, key));
        sqrtprec_vec = @(u_vec) vec(sqrtprec(mat(u_vec)));

        eigs_opts.issym = vec_mat_are_isometries;
        eigs_opts.isreal = true;
        lambdas = eigs(@(u_vec) ...
                      sqrtprec_vec(hess_vec(sqrtprec_vec(u_vec))), ...
                      n, dim, 'LM', eigs_opts);

    elseif canGetPrecon(problem)

        % There is a preconditioner, but we don't have its square root:
        % deal with the non-symmetric composition Precon o Hessian.

        prec = @(u_mat) tgt(getPrecon(problem, x, tgt(u_mat), storedb, key));
        prec_vec = @(u_vec) vec(prec(mat(u_vec)));
        % prec_inv_vec = @(u_vec) pcg(prec_vec, u_vec);

        eigs_opts.issym = false;
        eigs_opts.isreal = true;
        lambdas = eigs(@(u_vec) prec_vec(hess_vec(u_vec)), ...
                       n, dim, 'LM', eigs_opts);

    else

        error('No preconditioner is available in the problem structure.');

    end

    lambdas = sort(lambdas);

end

author: NicolasBoumal

repo: manopt


text:
classdef TestFastNlMeansDenoisingMulti
    %TestFastNlMeansDenoisingMulti

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','lena.jpg');
    end

    methods (Static)
        function test_1
            img = cv.imread(TestFastNlMeansDenoisingMulti.im, ...
                'Grayscale',true, 'ReduceScale',2);
            imgs = repmat({img},1,5);
            out = cv.fastNlMeansDenoisingMulti(imgs, 2, 3, 'H',20);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_2
            % we use IMNOISE Image Processing Toolbox
            if ~mexopencv.require('images')
                error('mexopencv:testskip', 'toolbox');
            end

            img = cv.imread(TestFastNlMeansDenoisingMulti.im, ...
                'Grayscale',true, 'ReduceScale',2);
            imgs = cell(1,5);
            for i=1:numel(imgs)
                imgs{i} = imnoise(img, 'gaussian');
            end
            out = cv.fastNlMeansDenoisingMulti(imgs, 2, 3, 'H',20);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_3
            files = dir(fullfile(mexopencv.root(),'test','left0*.jpg'));
            imgs = cell(1, numel(files));
            for i=1:numel(files)
                imgs{i} = cv.imread(fullfile(mexopencv.root(),'test',files(i).name), ...
                    'Grayscale',true, 'ReduceScale',2);
            end
            out = cv.fastNlMeansDenoisingMulti(imgs, 6-1, 3);
            validateattributes(out, {class(imgs{1})}, {'size',size(imgs{1})});
        end

        function test_error_argnum
            try
                cv.fastNlMeansDenoisingMulti();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [OutputError,flir_vid] = PIVlab_capture_flir_synced_capture(flir_vid,nr_of_images,do_realtime,ROI_live,frame_nr_display)

hgui=getappdata(0,'hgui');
image_handle_flir=getappdata(hgui,'image_handle_flir');
OutputError=0;

flir_frames_to_capture = nr_of_images*2;

%% capture data

while flir_vid.FramesAcquired < (flir_frames_to_capture) &&  getappdata(hgui,'cancel_capture') ~=1
	ima = image_handle_flir.CData;
	set(frame_nr_display,'String',['Image nr.: ' int2str(round(flir_vid.FramesAcquired/2))]);
	drawnow limitrate
	%% sharpness indicator
	sharpness_enabled = getappdata(hgui,'sharpness_enabled');
	if sharpness_enabled == 1 % sharpness indicator
		textx=1240;
		texty=950;
		[~,~] = PIVlab_capture_sharpness_indicator (ima,textx,texty);
	else
		delete(findobj('tag','sharpness_display_text'));
	end
	crosshair_enabled = getappdata(hgui,'crosshair_enabled');
	if crosshair_enabled == 1 %cross-hair
		%% cross-hair
		locations=[0.15 0.5 0.85];
		half_thickness=1;
		brightness_incr=101;
		ima_ed=ima;
		old_max=max(ima(:));
		for loca=locations
			%vertical
			ima_ed(:,round(size(ima,2)*loca)-half_thickness:round(size(ima,2)*loca)+half_thickness)=ima_ed(:,round(size(ima,2)*loca)-half_thickness:round(size(ima,2)*loca)+half_thickness)+brightness_incr;
			%horizontal
			ima_ed(round(size(ima,1)*loca)-half_thickness:round(size(ima,1)*loca)+half_thickness,:)=ima_ed(round(size(ima,1)*loca)-half_thickness:round(size(ima,1)*loca)+half_thickness,:)+brightness_incr;
		end
		ima_ed(ima_ed>old_max)=old_max;
		set(image_handle_flir,'CData',ima_ed);
	end
	%% HISTOGRAM
	if getappdata(hgui,'hist_enabled')==1
		if isvalid(image_handle_flir)
			hist_fig=findobj('tag','hist_fig');
			if isempty(hist_fig)
				hist_fig=figure('numbertitle','off','MenuBar','none','DockControls','off','Name','Live histogram','Toolbar','none','tag','hist_fig','CloseRequestFcn', @HistWindow_CloseRequestFcn);
			end
			if ~exist ('old_hist_y_limits','var')
				old_hist_y_limits =[0 35000];
			else
				if isvalid(hist_obj)
					old_hist_y_limits=get(hist_obj.Parent,'YLim');
				end
			end
			hist_obj=histogram(ima(1:2:end,1:2:end),'Parent',hist_fig,'binlimits',[0 255]);
		end
		%lowpass hist y limits for better visibility
		if ~exist ('new_hist_y_limits','var')
			new_hist_y_limits =[0 35000];
		end
		new_hist_y_limits=get(hist_obj.Parent,'YLim');
		set(hist_obj.Parent,'YLim',(new_hist_y_limits*0.5 + old_hist_y_limits*0.5))
	else
		hist_fig=findobj('tag','hist_fig');
		if ~isempty(hist_fig)
			close(hist_fig)
		end
	end
	drawnow limitrate;
	%% Autofocus
	%% Lens control
	%to be implemented...
end

stoppreview(flir_vid)
stop(flir_vid);
set(frame_nr_display,'String',['Image nr.: ' int2str(round(flir_vid.FramesAcquired/2))]);
drawnow;

function HistWindow_CloseRequestFcn(hObject,~)
hgui=getappdata(0,'hgui');
setappdata(hgui,'hist_enabled',0);
try
	delete(hObject);
catch
	delete(gcf);
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function f = flops_solve_tri(n,m,k)
% FLOPS_SOLVE_TRI   Flops for triangular left division.
% FLOPS_SOLVE_TRI(T,b) returns the number of flops for solve_tri(T,b).
% FLOPS_SOLVE_TRI(n,m,k) returns the number of flops for
% solve_tril(tril(rand(n,m)),rand(n,k)).
%
% Example:  (n=2,m=2,k=1)
%  [g;h] = [a 0; b c]\[e; f]  has
%  g = e/a
%  h = (f - b*g)/c
%  which is 2 multiply+add and 2 divisions = 18 flops.

if nargin == 2
	T = n;
	b = m;
  f = flops_solve_tri(rows(T),cols(T),cols(b));
  return;
end
if n ~= m
	error('n ~= m case is not implemented');
end
% lower triangular case:
% number of multiplies+adds is
% sum(i=1..n) sum(k=1..i-1) 2 = sum(i=1..n) 2*(i-1) = n^2-n
% number of divides is n
f = (n*n + n*(flops_div-1))*k;

author: andrewssobral

repo: lrslibrary


text:
function [EI,Pth,Qthth] = lme_EI(X,Zcols,W,CBhat,SIGMA,L,phi,ni)
% [EI,Pth,Qthth] = lme_EI(X,Zcols,W,CBhat,SIGMA,L,phi,ni)
%
% Expected information matrix for the restricted log-likelihood.
%
% Input
% X: Ordered design Matrix (according to time for each subject).
% Zcols: Vector with the indices of the colums of X that will be considered
% as random effects.
% W: Inverses of the estimated marginal covariance matrices for each
% subject stacked in W.
% CBhat: Asymptotic covariance matrix of the fixed effects.
% SIGMA: Estimated marginal covariance matrices for each subject
% stacked in SIGMA.
% L: Cholesky factor of the covariance matrix of the random effects (D).
% phi: Within-subject standard deviation of the errors.
% ni: Vector whose entries are the number of repeated measures for each
% subject in the study (ordered according to X).
%
% Output
% EI: Expected information matrix.
% Pth,Qthth: Matrices that are useful for inferences on the fixed effects.
%
% $Revision: 1.2 $  $Date: 2015/01/06 17:14:57 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:14:57 $
%    $Revision: 1.2 $
% Reference: Kenward MG and Roger JH, 1997. Small sample inference for fixed
% effects from restricted maximum likelihood. Biometrics,Vol. 53, No.3.
%
m = length(ni);
n = sum(ni);
q = length(Zcols);
nth = q*(q+1)/2+1;
p = size(X,2);
Pth = zeros(nth,p,p);
Qthth = zeros(nth,nth,p,p);
Der = zeros(nth,n,size(W,2));
%Computation of the first order derivatives of W
jk = 0;
for k=1:q
    for j=1:k
        jk = jk + 1;
        posi = 1;
        for i=1:m
            posf = posi+ni(i)-1;
            Zi = X(posi:posf, Zcols);
            Zki = Zi(:,k);
            Mjki = Zki*L(j,:)*Zi';
            Mjki = Mjki + Mjki';
            Wi = W(posi:posf,1:ni(i));
            Der(jk,posi:posf,1:ni(i)) = -Wi*Mjki*Wi;
            posi = posf+1;
        end;
    end;
end;
posi = 1;
for i=1:m
    posf = posi+ni(i)-1;
    Wi = W(posi:posf,1:ni(i));
    Der(nth,posi:posf,1:ni(i)) = -2*phi*Wi*Wi;
    posi = posf+1;
end;
%Computation of Pis,Qijs and the expected information matrix EI.
for j=1:nth
    posi = 1; Pj = 0;
    Bj = squeeze(Der(j,:,:));
    for i=1:m
        posf = posi+ni(i)-1;
        Pj = Pj + X(posi:posf,:)'*Bj(posi:posf,1:ni(i))*X(posi:posf,:);
        posi = posf+1;
    end;
    Pth(j,:,:) = Pj;
end;
EI = zeros(nth,nth);
for k=1:nth
    Bk = squeeze(Der(k,:,:));
    Pk = squeeze(Pth(k,:,:));
    for j=1:k
        Bj = squeeze(Der(j,:,:));
        Pj = squeeze(Pth(j,:,:));
        posi = 1; Qkj = 0;
        traceBkj = 0;
        for i=1:m
            posf = posi+ni(i)-1;
            Vi = SIGMA(posi:posf,1:ni(i));
            Bkji = Bk(posi:posf,1:ni(i))*Vi*Bj(posi:posf,1:ni(i));
            traceBkj = traceBkj + trace(Bkji*Vi);
            Qkji = X(posi:posf,:)'*Bkji*X(posi:posf,:);
            Qkj = Qkj + Qkji;
            posi = posf+1;
        end;
        Qthth(k,j,:,:) = Qkj;
        Qthth(j,k,:,:) = Qkj;
        EI(k,j) = 0.5*(traceBkj - trace(CBhat*(2*Qkj-Pk*CBhat*Pj)));
        EI(j,k) = EI(k,j);
    end;
end;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function gpmf = gpmf_constant(varargin)
%GPMF_CONSTANT  Create a constant mean function
%
%  Description
%    GPMF = GPMF_CONSTANT('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates constant mean function structure in which the named
%    parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    GPMF = GPMF_CONSTANT(GPMF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a mean function structure with the named parameters
%    altered with the specified values.
%
%    Parameters for constant mean function
%      constant          - constant value for the constant
%                          base function (default 1)
%      prior_mean        - prior mean (scalar or vector) for base
%                          functions' weight prior (default 0)
%      prior_cov         - prior covariances (scalar or vector)
%                          for base functions' prior corresponding
%                          each selected input dimension. In
%                          multiple dimension case prior_cov is a
%                          struct containing scalars or vectors.
%                          The covariances must all be either
%                          scalars (diagonal cov.matrix) or
%                          vectors (for non-diagonal cov.matrix)
%                          (default 100)
%
%  See also
%    GP_SET, GPMF_LINEAR, GPMF_SQUARED
%
% Copyright (c) 2010 Tuomas Nikoskinen
% Copyright (c) 2011 Jarno Vanhatalo

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

  ip=inputParser;
  ip.FunctionName = 'GPMF_CONSTANT';
  ip.addOptional('gpmf', [], @isstruct);
  ip.addParamValue('constant',1, @(x) isvector(x));
  ip.addParamValue('prior_mean',0, @(x) isvector(x));
  ip.addParamValue('prior_cov',100, @(x) isvector(x));
  ip.addParamValue('mean_prior', [], @isstruct);
  ip.addParamValue('cov_prior', [], @isstruct);
  ip.parse(varargin{:});
  gpmf=ip.Results.gpmf;

  if isempty(gpmf)
    % Initialize a mean function
    init=true;
    gpmf.type = 'gpmf_constant';
  else
    % Modify a mean function
    if ~isfield(gpmf,'type') && isequal(gpmf.type,'gpmf_constant')
      error('First argument does not seem to be a constant mean function')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('type',ip.UsingDefaults)
    gpmf.constant = ip.Results.constant;
  end
  if init || ~ismember('prior_mean',ip.UsingDefaults)
    gpmf.b=ip.Results.prior_mean(:)';
  end
  if init || ~ismember('prior_mean',ip.UsingDefaults)
    gpmf.B=ip.Results.prior_cov(:)';
  end
  if init || ~ismember('mean_prior',ip.UsingDefaults)
    gpmf.p.b=ip.Results.mean_prior;
  end
  if init || ~ismember('cov_prior',ip.UsingDefaults)
    gpmf.p.B=ip.Results.cov_prior;
  end
  if init
    % Set the function handles to the nested functions
    gpmf.fh.geth = @gpmf_geth;
    gpmf.fh.pak = @gpmf_pak;
    gpmf.fh.unpak = @gpmf_unpak;
    gpmf.fh.lp = @gpmf_lp;
    gpmf.fh.lpg = @gpmf_lpg;
    gpmf.fh.recappend = @gpmf_recappend;
  end

end

function h = gpmf_geth(gpmf, x)
%GPMF_GETH  Calculate the base function values for a given input.
%
%  Description
%    H = GPMF_GETH(GPMF,X) takes in a mean function structure
%    GPMF and inputs X. The function returns a row vector of
%    length(X) containing the constant value which is by default
%    1.

  constant=gpmf.constant;
  if isscalar(constant)
    h = repmat(constant,1,length(x(:,1)));
  else
    h = constant(:)';
  end

end

function [w, s, h] = gpmf_pak(gpmf, w)
%GPMF_PAK  Combine GP mean function parameters into one vector
%
%  Description
%    W = GPMF_PAK(GPMF) takes a mean function
%    structure GPMF and combines the mean function
%    parameters and their hyperparameters into a single row
%    vector W.
%
%       w = [ log(gpmf.b)
%             (hyperparameters of gpmf.b)
%             log(gpmf.B)
%             (hyperparameters of gpmf.B)]'
%
%  See also
%    GPMF_UNPAK

  w = []; s = {}; h=[];
  if ~isempty(gpmf.p.b)
    w = gpmf.b;
    if numel(gpmf.b)>1
      s = [s; sprintf('gpmf_constant.b x %d',numel(gpmf.b))];
    else
      s = [s; 'gpmf_constant.b'];
    end
    h = [h -1.*ones(1,numel(gpmf.b))];
    % Hyperparameters of b
    [wh, sh, hh] = gpmf.p.b.fh.pak(gpmf.p.b);
    w = [w wh];
    s = [s; sh];
    h = [h -1-hh];
  end

  if ~isempty(gpmf.p.B)
    w = [w log(gpmf.B)];
    if numel(gpmf.B)>1
      s = [s; sprintf('log(gpmf_constant.B x %d)',numel(gpmf.B))];
    else
      s = [s; 'log(gpmf_constant.B)'];
    end
    h = [h -1.*ones(1,numel(gpmf.B))];
    % Hyperparameters of B
    [wh, sh, hh] = gpmf.p.B.fh.pak(gpmf.p.B);
    w = [w wh];
    s = [s; sh];
    h = [h -1-hh];
  end

end

function [gpmf, w] = gpmf_unpak(gpmf, w)
%GPMF_UNPAK  Sets the mean function parameters into the structure
%
%  Description
%    [GPMF, W] = GPMF_UNPAK(GPMF, W) takes a covariance
%    function structure GPMF and a hyper-parameter vector W, and
%    returns a mean function structure identical to the
%    input, except that the covariance hyper-parameters have been
%    set to the values in W. Deletes the values set to GPMF from
%    W and returns the modified W.
%
%    Assignment is inverse of
%       w = [ log(gpmf.b)
%             (hyperparameters of gpmf.b)
%             log(gpmf.B)
%             (hyperparameters of gpmf.B)]'
%
%  See also
%   GPMF_PAK

  gpp=gpmf.p;

  if ~isempty(gpp.b)
    i2=length(gpmf.b);
    i1=1;
    gpmf.b = w(i1:i2);
    w = w(i2+1:end);

    % Hyperparameters of b
    [p, w] = gpmf.p.b.fh.unpak(gpmf.p.b, w);
    gpmf.p.b = p;
  end

  if ~isempty(gpp.B)
    i2=length(gpmf.B);
    i1=1;
    gpmf.B = exp(w(i1:i2));
    w = w(i2+1:end);

    % Hyperparameters of B
    [p, w] = gpmf.p.B.fh.unpak(gpmf.p.B, w);
    gpmf.p.B = p;
  end

end

function lp = gpmf_lp(gpmf)
%GPMF_SEXP_LP  Evaluate the log prior of mean function parameters
%
%  Description
%
%  See also

% Evaluate the prior contribution to the error. The parameters that
% are sampled are transformed, e.g., W = log(w) where w is all
% the "real" samples. On the other hand errors are evaluated in
% the W-space so we need take into account also the Jacobian of
% transformation, e.g., W -> w = exp(W). See Gelman et al. (2013),
% Bayesian Data Analysis, third edition, p. 21.
  lp = 0;
  gpp=gpmf.p;

  if ~isempty(gpmf.p.b)
    lp = lp + gpp.b.fh.lp(gpmf.b, ...
                   gpp.b);
  end

  if ~isempty(gpp.B)
    lp = lp + gpp.B.fh.lp(gpmf.B, ...
                   gpp.B) +sum(log(gpmf.B));
  end
end

function [lpg_b, lpg_B] = gpmf_lpg(gpmf)
%GPMF_SEXP_LPG  Evaluate gradient of the log prior with respect
%               to the parameters.
%
%  Description
%    LPG = GPMF_SEXP_LPG(GPMF) takes a mean function
%    structure GPMF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters.
%
%  See also
%    GPMF_SEXP_PAK, GPMF_SEXP_UNPAK, GPMF_SEXP_LP, GP_G

  lpg_b=[]; lpg_B=[];
  gpp=gpmf.p;

  if ~isempty(gpmf.p.b)
    lll = length(gpmf.b);
    lpgs = gpp.b.fh.lpg(gpmf.b, gpp.b);
    lpg_b = [lpgs(1:lll) lpgs(lll+1:end)];  %.*gpmf.b+1
  end

  if ~isempty(gpmf.p.B)
    lll = length(gpmf.B);
    lpgs = gpp.B.fh.lpg(gpmf.B, gpp.B);
    lpg_B = [lpgs(1:lll).*gpmf.B+1 lpgs(lll+1:end)];
  end
end

function recmf = gpmf_recappend(recmf, ri, gpmf)
%RECAPPEND  Record append
%
%  Description
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

% Initialize record
  if nargin == 2
    recmf.type = 'gpmf_constant';

    % Initialize parameters
    recmf.b= [];
    recmf.B = [];

    % Set the function handles
    recmf.fh.geth = @gpmf_geth;
    recmf.fh.pak = @gpmf_pak;
    recmf.fh.unpak = @gpmf_unpak;
    recmf.fh.lp = @gpmf_lp;
    recmf.fh.lpg = @gpmf_lpg;
    recmf.fh.recappend = @gpmf_recappend;

    recmf.p=[];
    recmf.p.b=[];
    recmf.p.B=[];
    if isfield(ri.p,'b') && ~isempty(ri.p.b)
      recmf.p.b = ri.p.b;
    end
    if ~isempty(ri.p.B)
      recmf.p.B = ri.p.B;
    end
    return
  end

  gpp = gpmf.p;

  % record magnSigma2
  if ~isempty(gpmf.b)
    recmf.b(ri,:)=gpmf.b;
    if ~isempty(recmf.p.b)
      recmf.p.b = gpp.b.fh.recappend(recmf.p.b, ri, gpmf.p.b);
    end
  elseif ri==1
    recmf.b=[];
  end

  if ~isempty(gpmf.B)
    recmf.B(ri,:)=gpmf.B;
    if ~isempty(recmf.p.B)
      recmf.p.B = gpp.B.fh.recappend(recmf.p.B, ri, gpmf.p.B);
    end
  elseif ri==1
    recmf.B=[];
  end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function obj = modelBatchCopyFiles(obj)
	% Copies files from folders to a specific other folder.
	% Biafra Ahanonu
	% 2015.09.01 [21:36:38]
	%
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

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

	% get user input
	movieSettings = inputdlg({...
			'start:end frames (leave blank for all)',...
			'file regexp:',...
			'copy to specific folder (leave blank to copy to same folder)',...
			'input HDF5 dataset name',...
			'output HDF5 dataset name',...
			'analyzing movie files (1 = yes, 0 = no):',...
			'back-up directory name (if copying into same folder)',...
			'name of new extension for movie (include leading dot, leave blank if keep the same)',...
			'append folder base string (blank if no)'...
		},...
		'copy files to /archive/ folder',1,...
		{...
			'1:4000',...
			obj.fileFilterRegexp,...
			'',...
			obj.inputDatasetName,...
			obj.inputDatasetName,...
			'1',...
			'archive',...
			'',...
			''...
		}...
	);
	setNo = 1;
	frameList = str2num(movieSettings{setNo});setNo = setNo+1;
	fileFilterRegexp = movieSettings{setNo};setNo = setNo+1;
	saveSpecificFolder = movieSettings{setNo};setNo = setNo+1;
	inputDatasetName =  movieSettings{setNo};setNo = setNo+1;
	outputDatasetName =  movieSettings{setNo};setNo = setNo+1;
	analyzeMovieFiles =  str2num(movieSettings{setNo});setNo = setNo+1;
	backupDirName =  movieSettings{setNo};setNo = setNo+1;
	newExtensionName =  movieSettings{setNo};setNo = setNo+1;
	useBaseFolderName =  movieSettings{setNo};setNo = setNo+1;

	[fileIdxArray idNumIdxArray nFilesToAnalyze nFiles] = obj.getAnalysisSubsetsToAnalyze();

	% loop over all directories and copy over a specified chunk of the movie file
	for thisFileNumIdx = 1:nFilesToAnalyze
			thisFileNum = fileIdxArray(thisFileNumIdx);
			obj.fileNum = thisFileNum;
			display(repmat('=',1,21))
			display([num2str(thisFileNumIdx) '/' num2str(nFilesToAnalyze) ' (' num2str(thisFileNum) '/' num2str(nFiles) '): ' obj.fileIDNameArray{obj.fileNum}]);
			% =====================
			switch analyzeMovieFiles
				case 1
					% load the movie
					movieList = getFileList(obj.inputFolders{obj.fileNum}, fileFilterRegexp);
					% movieList = movieList{1};
					nMovies = length(movieList);
					for movieNo = 1:nMovies
						moviePath = movieList{movieNo};
						if isempty(frameList)
							frameListTmp = frameList;
						else
							movieDims = loadMovieList(moviePath,'convertToDouble',0,'frameList',[],'inputDatasetName',inputDatasetName,'treatMoviesAsContinuous',1,'loadSpecificImgClass','single','getMovieDims',1);
							frameListTmp = frameList;
							% remove frames that are too large
							frameListTmp(frameListTmp>=movieDims.z) = [];
						end

						% create directory in other root
						if isempty(saveSpecificFolder)
							% move movie to archive folder if copying to the same folder then move file
							archivePath = [obj.inputFolders{obj.fileNum} filesep backupDirName];
							mkdir(archivePath);
							[PATHSTR,NAME,EXT] = fileparts(moviePath);
							archivePathFile = [archivePath filesep NAME EXT];
							display(['moving' num2str(movieNo) '/' num2str(nMovies) ': ' moviePath ' TO ' archivePathFile]);
							% movefile(moviePath,archivePathFile);
							if ispc
								dos(['move ' moviePath ' ' archivePathFile]);
							elseif isunix
								unix(['mv ' moviePath ' ' archivePathFile]);
							end
							newPathFile = moviePath;
							moviePath = archivePathFile;
						else
							[PATHSTR,NAME,EXT] = fileparts(obj.inputFolders{obj.fileNum});
							newPath = [saveSpecificFolder filesep NAME];
							mkdir(newPath)
							[PATHSTR,NAME,EXT] = fileparts(moviePath);
							newPathFile = [newPath filesep NAME EXT];
						end

						[primaryMovie] = loadMovieList(moviePath,'convertToDouble',0,'frameList',frameListTmp(:),'inputDatasetName',inputDatasetName);

						if ~isempty(newExtensionName)
							[PATHSTR,NAME,EXT] = fileparts(newPathFile);
							newPathFile = [PATHSTR filesep NAME newExtensionName];
						end

						if ~isempty(useBaseFolderName)
							[PATHSTR,NAME,EXT] = fileparts(newPathFile);
							newPathFile = [PATHSTR filesep obj.folderBaseSaveStr{obj.fileNum} '_' NAME EXT];
						end

						% save the file in the new location
						[output] = writeHDF5Data(primaryMovie,newPathFile,'datasetname',outputDatasetName);
					end
					% body
				case 0
					% load the movie
					fileList = getFileList(obj.inputFolders{obj.fileNum}, fileFilterRegexp);
					nFiles = length(fileList);
					for fileNo = 1:nFiles
						filePath = fileList{fileNo};
						% create directory in other root
						if isempty(saveSpecificFolder)
							% move files to folder in current directory
							archivePath = [obj.inputFolders{obj.fileNum} filesep backupDirName];
							mkdir(archivePath);
							[PATHSTR,NAME,EXT] = fileparts(filePath);
							archivePathFile = [archivePath filesep NAME EXT];
							display(['moving' num2str(fileNo) '/' num2str(nFiles) ': ' filePath ' TO ' archivePathFile]);
							% movefile(filePath,archivePathFile);
							if ispc
								dos(['move ' filePath ' ' archivePathFile]);
							elseif isunix
								unix(['mv ' filePath ' ' archivePathFile]);
							end
							newPathFile = filePath;
							filePath = archivePathFile;
						else
							[PATHSTR,NAME,EXT] = fileparts(obj.inputFolders{obj.fileNum});
							newPath = [saveSpecificFolder filesep NAME];
							mkdir(newPath)
							[PATHSTR,NAME,EXT] = fileparts(filePath);
							newPathFile = [newPath filesep NAME EXT];
							display(['copying' num2str(fileNo) '/' num2str(nFiles) ': ' filePath ' TO ' newPathFile]);
							if ispc
								dos(['copy ' filePath ' ' newPathFile]);
							elseif isunix
								unix(['cp ' filePath ' ' newPathFile]);
							end
						end
					end
				otherwise
					% body
			end
	end

end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [marginal, pot] = marginal_family(engine, query)
% MARGINAL_NODES Compute the marginal on the family of the specified query node (belprop)
% [marginal, pot] = marginal_family(engine, query)

pot = engine.marginal_domains{query};
marginal = pot_to_marginal(pot);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [subplot_handles,LabelsFontSize] = subplotplus(subplot_array,varargin)
% SUBPLOTPLUS Enhanced subplots creation function.
%
%   [subplot_handles,LabelsFontSize] = subplotplus(subplot_array,varargin)
%   is an enhanced alternative to the inherit subplot() function of
%   MATLAB. Using subplotplus() almost any layout of subplots can be
%   created in a particular figure. The desired subplot layout is described
%   using a complex cell array where each "atomic" cell, within the cell
%   array corresponds to a subplot. The function scans the content of the
%   input cell array and build the desired matching subplot layout in the
%   figure. The size, alignment and font size of the subplots is
%   automatically adjusted for best appearance. The function is recursive
%   in nature so any level of subplot divisions can be created. The
%   returned arguments are the created subplot handles as well as the
%   optimal label (X and Y) font size that should be used for a specific
%   subplot for best appearance.
%
%   Input arguments:
%   + subplot_array: A cell array that describes the desired subplot
%   layput. For example, a 1x2 subplot layout where the right-end subplot is
%   further divided down to a column of 3 additional subplots is described
%   by the following cell expression: {{[]},{{[]};{[]};{[]}}}. A group of
%   subplot cells, sharing the same "parent" cell, can be glued together in order
%   to save figure space using the '-g' option at the cell's content (instead of
%   empty group []). For example, if the 3 columned subplots of the
%   previous example are to be glued together, sharing the same X-axis, the
%   follwing cell expression should be used:
%   {{[]},{{['...-g...']};{['...-g...']};{['...-g...']}}}.
%   + varargin: only used for function recursion and should not be used by the user.
%
%   Output arguments:
%   + subplot_handles: A list of the created subplot handles by their order
%   of appearnace in the figure, starting from the left-top corner in a
%   clockwise manner.
%   + LabelsFontSize: A matching list of X/Y-axis optimal labels size for further
%   user work on the subplots.
%
%   Alon Geva
%   $Revision: 1.00 $  $Date: 12/01/2012 18:52:53 $

arraydim =  cell_dimension(subplot_array);

if (nargin == 1)
    clf;
    for i=1:arraydim,
        subplot(prod(arraydim),1,i);
    end
    subplot_handles = get(gcf,'children');
    subplot_location = [0.015 0 1-0.015 1]; %full plot region
    subplot_handles = flipud(subplot_handles);
    original_subplot_handles = subplot_handles;
    LabelsFontSize = zeros(length(original_subplot_handles),2);
else
    subplot_handles = varargin{1};
    subplot_location = varargin{2};
    LabelsFontSize = varargin{3};
    original_subplot_handles = -1;
end

subplotarraysize=size(subplot_array);
%i=1;
glue_FLAG = 0;
if ~iscell(subplot_array{1,1})
    if ischar(subplot_array{1,1})
        glue_FLAG = ~isempty(regexp(subplot_array{1,1},'-g','start'));
    end
elseif ~iscell(subplot_array{1,1}{1,1})
    if ischar(subplot_array{1,1}{1,1})
        glue_FLAG = ~isempty(regexp(subplot_array{1,1}{1,1},'-g','start'));
    end
end

if (glue_FLAG == 1)%(strcmp(subplot_array{1,1},'-g')==1)%(~isempty(regexp(subplot_array{1,1},'-g','start')))
    subplot_location(3) = subplot_location(3)*0.8;
    subplot_location(4) = subplot_location(4)*0.8;
    subplot_location(1) = subplot_location(1)+subplot_location(3)*0.1;
    subplot_location(2) = subplot_location(2)+subplot_location(4)*0.1;
end

for m=1:subplotarraysize(1),
    for n=1:subplotarraysize(2),
        %s = size(subplot_array{m,n});

        if (iscell(subplot_array{m,n}))
            recursion_FLAG = iscell(subplot_array{m,n}{1,1});
        else
            recursion_FLAG = 0;
        end

        if (~recursion_FLAG)%iscell(subplot_array{m,n}{1,1})%(s(1)<=1 && s(2)<=1)%
            M = subplotarraysize(1);
            N = subplotarraysize(2);
            if (~glue_FLAG)%(strcmp(subplot_array{1,1},'-g')~=1)%(isempty(regexp(subplot_array{1,1},'-g','start')))%
                set(subplot_handles(1),'position',[subplot_location(1)+(2*n-1)*subplot_location(3)/(2*N)-0.5*subplot_location(3)/N/(0.5/0.3347),...
                                                    subplot_location(2)+subplot_location(4)-(2*m-1)*subplot_location(4)/(2*M)-0.5*subplot_location(4)/M/(0.5/0.3412),...
                                                    subplot_location(3)/N/(0.5/0.3347),...
                                                    subplot_location(4)/M/(0.5/0.3412)]);
                Xfontfactor = subplot_location(3)/N/(0.5/0.3347);
                Yfontfactor = subplot_location(4)/M/(0.5/0.3412);
            else
                set(subplot_handles(1),'position',[subplot_location(1)+(2*n-1)*subplot_location(3)/(2*N)-0.5*subplot_location(3)/N/(0.5/0.5),...
                                                    subplot_location(2)+subplot_location(4)-(2*m-1)*subplot_location(4)/(2*M)-0.5*subplot_location(4)/M/(0.5/0.5),...
                                                    subplot_location(3)/N/(0.5/0.5),...
                                                    subplot_location(4)/M/(0.5/0.5)]);
                Xfontfactor = subplot_location(3)/N/(0.5/0.5);
                Yfontfactor = subplot_location(4)/M/(0.5/0.5);

                if (m == 1)
                    %set(subplot_handles(1),'XAxis','on');
                    set(subplot_handles(1),'XAxisLocation','top');
                elseif (m == subplotarraysize(1))
                    %set(subplot_handles(1),'XAxis','on');
                    set(subplot_handles(1),'XAxisLocation','bottom');
                else
                    set(subplot_handles(1),'XAxisLocation','bottom');
                end

                if (n == 1)
                    set(subplot_handles(1),'YAxisLocation','left');
                elseif (n == subplotarraysize(2))
                    set(subplot_handles(1),'YAxisLocation','right');
                else
                    set(subplot_handles(1),'YAxisLocation','right');
                end

%                 if (subplotarraysize(1) > 1)
%                     ylimmode = get(subplot_handles(1),'YLimMode');
%                     if (strcmp(ylimmode,'auto')==1)
%                         %ytick = get(subplot_handles(1),'YTick');
%                         ylim = get(subplot_handles(1),'YLim');
%                         set(subplot_handles(1),'YLimMode','manual','YLim',[ylim(1)-0.1*diff(ylim) ylim(2)+0.1*diff(ylim)]);
%                         %set(subplot_handles(1),'YTickMode','manual','YTick',ytick);
%                     end
%                 end
%
%                 if (subplotarraysize(2) > 1)
%                     xlimmode = get(subplot_handles(1),'XLimMode');
%                     if (strcmp(xlimmode,'auto')==1)
%                         %ytick = get(subplot_handles(1),'YTick');
%                         xlim = get(subplot_handles(1),'XLim');
%                         set(subplot_handles(1),'XLimMode','manual','XLim',[xlim(1)-0.1*diff(xlim) xlim(2)+0.1*diff(xlim)]);
%                         %set(subplot_handles(1),'YTickMode','manual','YTick',ytick);
%                     end
%                 end

            end

            set(subplot_handles(1),'XGrid','on','YGrid','on','Box','on','XMinorTick','on','YMinorTick','on','XMinorGrid','off','YMinorGrid','off');

            fontsize = 10+3*log10(min(Xfontfactor,Yfontfactor));
            if (fontsize > 10)
                fontsize = 10;
            elseif (fontsize < 1)
                fontsize = 1;
            end

            LabelsFontSize_tmp = 10+3*log10([Xfontfactor,Yfontfactor]);
            LabelsFontSize_tmp(LabelsFontSize_tmp > 10) = 10;
            LabelsFontSize_tmp(LabelsFontSize_tmp < 1) = 1;
            LabelsFontSize(length(LabelsFontSize)-length(subplot_handles)+1,:) = LabelsFontSize_tmp;

            set(subplot_handles(1),'FontSize',fontsize);
            subplot_handles = subplot_handles(2:end);
            %i = i + 1; %advancing subplots handle index
        else
            [subplot_handles,LabelsFontSize] = subplotplus(subplot_array{m,n},...
                                                        subplot_handles(1:end),...
                                                        [subplot_location(1)+(n-1)*subplot_location(3)/subplotarraysize(2),...
                                                         subplot_location(2)+subplot_location(4)-m*subplot_location(4)/subplotarraysize(1),...
                                                         subplot_location(3)/subplotarraysize(2),...
                                                         subplot_location(4)/subplotarraysize(1)],...
                                                         LabelsFontSize);
        end
    end
end

if (original_subplot_handles ~= -1) %root of recurrence
    subplot_handles = original_subplot_handles;
end

%%
function dimension = cell_dimension(C)

cellsize=size(C);
%p=prod(size(cell));
dimension=0;

for m=[1:cellsize(1)],
    for n=[1:cellsize(2)],
        s = size(C{m,n});
        if (s(1)<=1 && s(2)<=1)
            dimension = dimension + 1;
        else
            dimension = dimension + cell_dimension(C{m,n});
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function v = volume(SO3F,center,radius,varargin)
% ratio of orientations with a certain orientation
%
% Description
% The function 'volume' returns the ratio of an orientation that is close
% to an orientation (center) by a misorientation tolerance (radius) to the
% volume of the entire odf.
%
% Syntax
%   v = volume(odf,center,radius)
%   v = volume(odf,fibre,radius) % gives the volume with a fibre
%
% Input
%  odf    - @SO3Fun
%  center - @orientation
%  fibre  - @fibre
%  radius - double
%
% Options
%  resolution - resolution of discretization
%
% See also
% SO3Fun/fibreVolume SO3Fun/entropy SO3Fun/textureindex

if isa(center,'fibre')

  v = fibreVolume(SO3F,center.h,center.r,radius,varargin{:});

else

  % get resolution
  res = get_option(varargin,'RESOLUTION',min(1.25*degree,radius/30),'double');

  % discretisation
  if nargin > 3 && isa(varargin{1},'orientation')
    S3G = varargin{1};
  else
    S3G = equispacedSO3Grid(SO3F.CS,SO3F.SS,...
      'maxAngle',radius,'center',center,'resolution',res,varargin{:});
  end

  % estimate volume portion of odf space
  reference = 9897129 * 96 / numProper(SO3F.CS) / numProper(SO3F.SS);
  f = min(1,length(S3G) * (res / 0.25 / degree)^3 / reference);

  % eval odf
  if f == 0
    v = 0;
  else
    v = mean(eval(SO3F,S3G)) * f;
  end
end

end

author: mtex-toolbox

repo: mtex

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


text:
function f = f01_f0 ( n, x, y )

%*****************************************************************************80
%
%% F01_F0 returns the value of function 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real F(N,1), the function values.
%
  f(1:n,1) = ...
      0.75 * exp ( - ( ( 9.0 * x(1:n,1) - 2.0 ).^2            ...
                     + ( 9.0 * y(1:n,1) - 2.0 ).^2 ) / 4.0 )  ...
    + 0.75 * exp ( - ( ( 9.0 * x(1:n,1) + 1.0 ).^2 ) / 49.0   ...
                     - ( 9.0 * y(1:n,1) + 1.0 ) / 10.0 )      ...
    + 0.5  * exp ( - ( ( 9.0 * x(1:n,1) - 7.0 ).^2            ...
                     + ( 9.0 * y(1:n,1) - 3.0 ).^2 ) / 4.0 )  ...
    - 0.2  * exp ( - (   9.0 * x(1:n,1) - 4.0 ).^2            ...
                     - ( 9.0 * y(1:n,1) - 7.0 ).^2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%RG_ADDTICKS Label spectral locus
%
% RG_ADDTICKS() adds wavelength ticks to the spectral locus.
%
% See also XYCOLOURSPACE.

% 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/>.

function rg_addticks(lam1, lam2, lamd)

    % well spaced points around the locus
    lambda = [460:10:540];
    lambda = [lambda 560:20:600];

    rgb = cmfrgb(lambda*1e-9);
    r = rgb(:,1)./sum(rgb')';
    g = rgb(:,2)./sum(rgb')';
    hold on
    plot(r,g, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 6)
    hold off

    for i=1:numcols(lambda)
        text(r(i), g(i), sprintf('  %d', lambda(i)));
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [AD,Rank] = UpdateDiversityArchive(AD,Population,ND,R,Z,Xre,sigma_niche,Problem)
% Update the Diversity Archive

%------------------------------- Copyright --------------------------------
% Copyright 2017-2018 Yiping Liu
% Please contact {yiping0liu@gmail.com} if you have any problem.
%--------------------------------------------------------------------------

    %% Non-dominated sorting
    Population = [Population,AD];
    [FrontNo,MaxFNo] = NDSort(Population.objs,ND);

    %% Select the solutions in the last front
    if length(find(FrontNo<=MaxFNo))== ND
       Next = FrontNo<=MaxFNo;
    else
       Next = FrontNo < MaxFNo ;
       NQ     = ND-sum(Next);
       Last   = find(FrontNo==MaxFNo);
       Choose = LastSelection(Population(Last).objs,Population(Last).decs,NQ,R,Z,Xre,sigma_niche,Problem); %Algorithm 4, lines 5-30
       Next(Last(Choose)) = true;
    end

    %% Population for next generation
    AD = Population(Next);
    Rank = FrontNo(Next);
end

function Choose = LastSelection(PopObj,PopDec,NQ,R,Z,Xre,sigma_niche,Problem)
% Select solutions with good diversity in the last front
    N      = size(PopObj,1);
    NR     = size(R,1);

    %% Normalize
    PopObj = PopObj - repmat(Z,N,1);
    PopDec = (PopDec - repmat(Problem.lower,N,1))./repmat(Problem.upper-Problem.lower,N,1);

    %% Calculate distance between every solution and reference point in the objective space
    theta = pdist2(R,PopObj,'cosine');

    %% Calculate distance between every two solutions in the reminder decsion subspace
    d = pdist2(PopDec(:,Xre),PopDec(:,Xre),'chebychev');

    %% Cluster
    [thmin,label] = min(theta,[],1);
    C1 = false(NR,N);
    C2 = false(NR,N);
    for j = 1:NR
        member = find(label==j);
        member1 = label==j;
        [~,temp] = sort(thmin(member));
        for i = temp
           if any(d(member(i),and(C1(j,:),member1))<sigma_niche)
               C2(j,member(i))=true;
           else
               C1(j,member(i))=true;
           end
        end
    end

    %% Make selected solution == NQ
    while sum(sum(C1))>NQ
        cmax = max(sum(C1,2));
        jmax = sum(C1,2)== cmax;
        temp1 = find(sum(C1(jmax,:),1)>0);
        [~,xmax] = max(thmin(temp1));
        xmax=temp1(xmax);
        C1(label(xmax),xmax) = false;
    end
    while sum(sum(C1))<NQ
        c2 = sum(C2,2)>0;
        cmin = min(sum(C1(c2,:),2));
        jmin = and(sum(C1,2)==cmin,c2);
        temp1 = find(sum(C2(jmin,:),1)>0);
        [~,xmin] = min(thmin(temp1));
        xmin=temp1(xmin);
        C2(label(xmin),xmin) = false;
        C1(label(xmin),xmin) = true;
    end
    Choose = sum(C1)>0;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [x,g,xn,gg] = kmeanhar(d,k,l,e,x0)
%KMEANS Vector quantisation using K-harmonic means algorithm [X,G,XN,GG]=(D,K,L,E,X0)
%
%  Inputs:
%
%    D(N,P)  contains N data vectors of dimension P
%    K       is number of centres required
%    L       integer portion is max loop count, fractional portion
%            gives stopping threshold as fractional reduction in performance criterion
%    E       is exponent in the cost function. Significantly faster if this is an even integer. [default 4]
%    X0(K,P) are the initial centres (optional)
%            Alternatively, X0 can be a character determining the initialization method:
%                'f'    Initialize with K randomly selected data points [default]
%                'p'    Initialize with centroids and variances of random partitions
%
%  Outputs:
%
%    X(K,P)  is output row vectors
%    G       is the final performance criterion value (normalized by N)
%    XN      nearest centre for each input point
%    GG(L+1) value of performance criterion before each iteration and at end
%
% The k-harmonic means algorithm selects K cluster centres to minimize
%                           sum_n(K/sum_k((d_n-x_k)^-e))
% where sum_n is over the N inputs points d_n and sum_k is over the K cluster centres x_k.
%
% It is often a good idea to scale the input data so that it has equal variance in each
% dimension before calling KMEANHAR so that approximately equal weight is given
% to each dimension in the distance calculation.

%  [1] Bin Zhang, "Generalized K-Harmonic Means - Boosting in Unsupervised Learning",
%      Hewlett-Packartd Labs, Technical Report HPL-2000-137, 2000 [Zhang2000]
%      http://www.hpl.hp.com/techreports/2000/HPL-2000-137.pdf

%  Bugs:
%      (1) Could use nested blocking to allow very large data arrays
%      (2) Could then allow incremental calling with partial data arrays (but messy)

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: kmeanhar.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% sort out the input arguments

if nargin<5
    x0='f';
    if nargin<4
        e=[];
        if nargin<3
            l=[];
        end
    end
end
if isempty(e)
    e=4;  % default value
end
if isempty(l)
    l=50+1e-3; % default value
end
sd=5;       % number of times we must be below threshold

% split into chunks if there are lots of data points

memsize=voicebox('memsize');
[n,p] = size(d);
nb=min(n,max(1,floor(memsize/(8*p*k))));    % block size for testing data points
nl=ceil(n/nb);                  % number of blocks

% initialize if X0 argument is not supplied

if ischar(x0)
    if k<n
        if any(x0=='p')                  % Initialize using a random partition
            ix=ceil(rand(1,n)*k);       % allocate to random clusters
            ix(rnsubset(k,n))=1:k;      % but force at least one point per cluster
            x=zeros(k,p);
            for i=1:k
                x(i,:)=mean(d(ix==i,:),1);
            end
        else                                % Forgy initialization: choose k random points [default]
            x=d(rnsubset(k,n),:);         % sample k centres without replacement
        end
    else
        x=d(mod((1:k)-1,n)+1,:);    % just include all points several times
    end
else
    x=x0;
end
eh=e/2;
th=l-floor(l);
l=floor(l)+(nargout>1);   % extra loop needed to calculate final performance value
if l<=0
    l=100;      % max number of iterations ever
end
if th==0
    th=-1;      % prevent any stopping if l has no fractional part
end
gg=zeros(l+1,1);
im=repmat(1:k,1,nb); im=im(:);

% index arrays for replication

wk=ones(k,1);
wp=ones(1,p);
% wn=ones(1,n);
%
% % Main calculation loop
%
% We have the following relationships to [1] where i and k index
% the data values and cluster centres respectively:
%
%   This program     [Zhang2000]                            Equation
%
%     d(i,:)            x_i                                 input data
%     x(k,:)            m_k                                 cluster centres
%     py(k,i)           (d_ik)^2
%     dm(i)'            d_i,min^2
%     pr(k,i)           (d_i,min/d_ik)^2
%     pe(k,i)           (d_i,min/d_ik)^p                    (7.6)
%     qik(k,i)          q_ik                                (7.2)
%     qk(k)             q_k                                 (7.3)
%     qik(k,i)./qk(k)   p_ik                                (7.4)
%     se(i)'            d_i,min^p * sumk(d_ik^-p)
%     xf(i)'            d_i,min^-2 / sumk(d_ik^-p)
%     xg(i)'            d_i,min^-(p+2) / sumk(d_ik^-p)^2

ss=sd+1;        % one extra loop at the start
g=0;                % dummy initial value of g
xn=zeros(n,1);
for j=1:l

    g1=g;                           % save old performance
    x1=x;                           % save old centres
    % first do partial chunk

    jx=n-(nl-1)*nb;
    ii=1:jx;
    kx=repmat(ii,k,1);
    km=repmat(1:k,1,jx);
    py=reshape(sum((d(kx(:),:)-x(km(:),:)).^2,2),k,jx);
    [dm,xn(ii)]=min(py,[],1);                 % min value in each column gives nearest centre
    dmk=dm(wk,:);                   % expand into a matrix
    dq=py>dmk;                      % update only these values
    pr=ones(k,jx);                   % leaving others at 1
    pr(dq)=dmk(dq)./py(dq);            % ratio of min(py)./py
    pe=pr.^eh;
    se=sum(pe,1);
    xf=dm.^(eh-1)./se;
    g=xf*dm.';                     % performance criterion (divided by k)
    xg=xf./se;
    qik=xg(wk,:).*pe.*pr;           % qik(k,i) is equal to q_ik in [Zhang2000]
    qk=sum(qik,2);
    xs=qik*d(ii,:);
    ix=jx+1;
    for il=2:nl
        jx=jx+nb;        % increment upper limit
        ii=ix:jx;
        kx=ii(wk,:);
        py=reshape(sum((d(kx(:),:)-x(im,:)).^2,2),k,nb);
        [dm,xn(ii)]=min(py,[],1);                 % min value in each column gives nearest centre
        dmk=dm(wk,:);                   % expand into a matrix
        dq=py>dmk;                      % update only these values
        pr=ones(k,nb);                   % leaving others at 1
        pr(dq)=dmk(dq)./py(dq);            % ratio of min(py)./py
        pe=pr.^eh;
        se=sum(pe,1);
        xf=dm.^(eh-1)./se;
        g=g+xf*dm.';                     % performance criterion (divided by k)
        xg=xf./se;
        qik=xg(wk,:).*pe.*pr;           % qik(k,i) is equal to q_ik in [Zhang2000]
        qk=qk+sum(qik,2);
        xs=xs+qik*d(ii,:);
        ix=jx+1;
    end
    gg(j)=g;
    x=xs./qk(:,wp);
    if g1-g<=th*g1
        ss=ss-1;
        if ~ss break; end  %  stop if improvement < threshold for sd consecutive iterations
    else
        ss=sd;
    end
end
gg=gg(1:j)*k/n;                       % scale and trim the performance criterion vector
g=g(end);
% gg' % *** DEBUIG ***
if nargout>1
    x=x1;                               % go back to the previous x values if G and/or XN value is output
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [poolfeat] = filelist_feature(dataset_name, filelist, feature, c)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

if(~isempty(dataset_name) && c.common_dictionary==0)
    c.cache=[c.cache '/' dataset_name '/'];
end

p = c.feature_config.(feature);
if(~isfield(p, 'dictionary'))
    c.feature_config.(feature).dictionary = build_dictionary(filelist, feature, c);
end

if(isempty(c.feature_config.(feature).dictionary))
    % Feature does not require LLC encoding + max pooling
    poolfeat = cell(length(filelist), 1);
    parfor i=1:length(filelist)
        img = imgread(filelist{i}, p);
        poolfeat{i} = extract_feature(feature, img, c);
    end
    poolfeat=cell2mat(poolfeat);
else
    % Feature requires LLC encoding + max pooling
    llcfeat = cell(length(filelist), 1);
    info = cell(length(filelist), 1);
    parfor i=1:length(filelist)
        img = imgread(filelist{i}, p);
        [llcfeat{i}, x, y, wid, hgt] = llc_feature(feature, img, c);
        info{i}.x = x; info{i}.y = y; info{i}.wid = wid; info{i}.hgt = hgt;
    end
    poolfeat = max_pooling(llcfeat, info, c.pool_region, p.pyramid_levels);
end

poolfeat = cast(poolfeat, c.precision);

author: adikhosla

repo: feature-extraction


text:
% Code used in Vidaurre et al. (2017) PNAS
%
% Detailed documentation and further examples can be found in:
% https://github.com/OHBA-analysis/HMM-MAR
% This pipeline must be adapted to your particular configuration of files.
%%%%%%%%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS AND FILE NAMES
mydir = '/home/diegov/MATLAB/';
addpath(genpath([ mydir 'HMM-MAR']))

K = 12; % no. states
repetitions = 3; % to run it multiple times (keeping all the results)
DirData = [mydir 'data/HCP/TimeSeries/hippocampus_data/'];
DirOut = [mydir 'HCP/HMM/out/hippocampus/'];
TR = 0.72;
use_stochastic = 1; % set to 1 if you have loads of data

N = 820; % no. subjects
f = cell(N,1); T = cell(N,1);
% each .mat file contains the data (ICA components) for a given subject,
% in a matrix X of dimension (4800time points by 50 ICA components).
% T{j} contains the lengths of each session (in time points)
for j=1:N
    f{j} = [DirData 's' num2str(j) '.mat'];
    load(f{j}); T{j} = [1200 1200 1200 1200];
end

options = struct();
options.K = K; % number of states
options.order = 0; % no autoregressive components
options.zeromean = 0; % model the mean
options.covtype = 'full'; % full covariance matrix
options.Fs = 1/TR;
options.verbose = 1;
options.standardise = 1;
options.inittype = 'HMM-MAR';
options.cyc = 500;
options.initcyc = 10;
options.initrep = 3;

% stochastic options
if use_stochastic
    options.BIGNbatch = round(N/30);
    options.BIGtol = 1e-7;
    options.BIGcyc = 500;
    options.BIGundertol_tostop = 5;
    options.BIGforgetrate = 0.7;
    options.BIGbase_weights = 0.9;
end

% We run the HMM multiple times
for r = 1:repetitions
    [hmm, Gamma, ~, vpath] = hmmmar(f,T,options);
    save([DirOut 'HMMrun_rep' num2str(r) '.mat'],'Gamma','vpath','hmm')
    disp(['RUN ' num2str(r)])
end

%% Pull out metastates

for r = 1:repetitions
    figure(r)
    load([DirOut 'HMMrun_rep' num2str(r) '.mat'],'Gamma','hmm')
    subplot(1,2,1) % Figure 2B
    GammaSessMean = squeeze(mean(reshape(Gamma,[1200 4 N K]),1));
    GammaSubMean = squeeze(mean(GammaSessMean,1));
    [~,pca1] = pca(GammaSubMean','NumComponents',1);
    [~,ord] = sort(pca1);
    imagesc(corr(GammaSubMean(:,ord))); colorbar
    subplot(1,2,2) % Figure 2A
    P = hmm.P;
    for j=1:K, P(j,j) = 0; P(j,:) = P(j,:) / sum(P(j,:));  end
    imagesc(P(ord,ord),[0 0.25]); colorbar
    axis square
    hold on
    for j=0:13
        plot([0 13] - 0.5,[j j] + 0.5,'k','LineWidth',2)
        plot([j j] + 0.5,[0 13] - 0.5,'k','LineWidth',2)
    end
    hold off
end

author: OHBA-analysis

repo: HMM-MAR


text:
function rotated_coords = rotatepoints(input_XY,center,anti_clockwise_angle, scale, varargin)
degree = 1; %Radians : degree = 0; Default is calculations in degrees

% Process the inputs
if length(varargin) ~= 0
    for n = 1:1:length(varargin)
        if strcmp(varargin{n},'degree')
            degree = 1;
        elseif strcmp(varargin{n},'radians')
            degree = 0;
        end
    end
    clear n;
end
[r,c] = size(input_XY);
if c ~= 2
    error('Not enough columns in coordinates XY ');
end
[r,c] = size(center);
if (r~=1 & c==2) | (r==1 & c~=2)
    error('Error in the size of the "center" matrix');
end

% Format the coordinate of the center of rotation
center_coord = input_XY;
center_coord(:,1) = center(1);
center_coord(:,2) = center(2);

% Turns the angles given to be such that the +ve is anti-clockwise and -ve is clockwise
anti_clockwise_angle = -1*anti_clockwise_angle;
% if in degrees, convert to radians because that's what the built-in functions use.
if degree == 1
    anti_clockwise_angle = deg2rad(anti_clockwise_angle);
end

%Produce the roation matrix
rotation_matrix = [cos(anti_clockwise_angle),-1*sin(anti_clockwise_angle);...
                   sin(anti_clockwise_angle),cos(anti_clockwise_angle)];
%Calculate the final coordinates
rotated_coords = scale*((input_XY-center_coord) * rotation_matrix) + center_coord;

end

author: jwyang

repo: face-alignment


text:
classdef ballfunv
%BALLFUNV   BALLFUNV class for representing vector-valued functions on the unit ball.
%
%  Class for approximating smooth vector-valued functions defined on the unit ball.
%
% BALLFUNV(F, G, H) constructs a BALLFUNV object representing the vector-valued
% function [F;G;H] on the unit ball. F, G, and H may be BALLFUN objects, function
% handles or scalars. If they are function handles, then those function handles
% should be vectorized.
%
% See also BALLFUN.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CLASS PROPERTIES:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties

        comp   % The three components Vx, Vy, Vz of a BALLFUNV

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )
        % Take the 3 components
        function F = ballfunv(varargin)

            % Return an empty BALLFUNV:
            if ( (nargin == 0) || isempty(varargin) )
                F.comp = {};
                return
            end

            % If the argument is a BALLFUNV, nothing to do:
            if ( isa(varargin{1}, 'ballfunv') )
                F = varargin{1};
                return
            end

            % BALLFUNV objects are vector-valued so complain if there
            % are less than 3 components:
            if ( numel(varargin) < 2 )
                error('BALLFUNV:ballfunv', ...
                    'Less than three components is not supported.')
            end

            % BALLFUNV objects cannot contain more than five components.
            % Complain if we have been given six or more.
            if ( numel(varargin) > 5 )
                error('BALLFUNV:ballfunv', ...
                    'More than four components is not supported.')
            end

            % Create a BALLFUNV from 3 ballfun or function handles
            if ( numel(varargin) == 3 )
                for jj = 1:3
                   if isa(varargin{jj}, 'ballfun') == 0
                       varargin{jj} = ballfun( varargin{jj} );
                   end
                end

                fh{1} = varargin{1};
                fh{2} = varargin{2};
                fh{3} = varargin{3};
            end

            % Create a BALLFUNV from 3 function handles in spherical
            % coordinates
            if ( numel(varargin) == 4 )
                if strcmp(varargin{4}, 'spherical') || strcmp(varargin{4}, 'polar')
                    fh{1} = ballfun(varargin{1},'spherical');
                    fh{2} = ballfun(varargin{2},'spherical');
                    fh{3} = ballfun(varargin{3},'spherical');
                else
                    error('BALLFUNV:ballfunv', ...
                        'Input arguments not supported')
                end
            end

            F.comp = fh;
        end
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% STATIC METHODS:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )
        varargout = PT2ballfunv(P,T);
    end

    methods ( Access = private, Static = true )

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function imOut = warpUV(imIn, u, v, x, y, interpolationMethod, template)
%
% function imOut = vWarpUV(imIn, u, v, x, y, interpolationMethod, template)
%
% DESCRIPTION
%   This function warps an input image according to the displacement
%     fields for all pixels (u and v) by using the MATLAB function interp2.
%     The interpolation method, e.g., bilinear interpolation or
%     nearest-neighbor interpolation, can be specified. Since this
%     function is called often, e.g., many times during the
%     registration process, the function tries to use stored
%     (persistent) matrices that are needed as arguments for
%     interp2 instead of recreating the matrices.
%
% INPUT
%   - imIn = image that is to be warped
%   - u = pixel displacement field in x direction
%   - v = pixel displacement field in y direction
%   - x is optional; can be given to the function to improve speed
%   - y is optional; can be given to the function to improve speed
%   - interpolationMetho is optional; 'nearest' is the default
%     'nearest':  edges better because of NaNs outside
%     'linear':   smoother inside
%   - template (optional): NaN regions in this image will be NaN in the
%     warped image as well; this can be the original phWedgeIm
%
% OUTPUT / RESULTS
%   - imOut = warped image
%
% AUTHOR:
%   Volker Maximillian Koch
%   vk@volker-koch.de
%
% DATE:
%   January - June 2001
%
% COMMENTS:
%   This function tries to avoid the calculation of x and y with meshgrid
%   since that takes some time. Therefore, these variables are stored
%   persistent.
%
% TO-DO:
%
% Update History:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent xStored;
persistent yStored;
persistent sizeStored;
if ~exist('x','var'),  x=[]; end
if ~exist('y','var'),  y=[]; end
if (isempty(x) | isempty(y))
    if isempty(sizeStored)
        sizeStored = [0 0];
    end
    if sizeStored == size(imIn)
        x = xStored;
        y = yStored;
    else
        [x,y] = meshgrid(1:size(imIn,2),1:size(imIn,1));
    end
end
xStored = x;
yStored = y;
sizeStored = size(imIn);
if eval('isempty(interpolationMethod)' , '1')
    interpolationMethod = '*linear';
end
imOut = zeros(size(imIn));
% we could use Bernd's version of interp2 (updateTinC).
%imOut = updateTinC(imIn, x+u, y+v);
imOut = interp2(imIn, x+u, y+v, interpolationMethod);
if eval('~isempty(template)' , '0')
    imOut(isnan(template)) = NaN;
end
return


text:
function [h] = hline(y, linetype)

if nargin < 2
    linetype={'r:'};
end

xlims = get(gca,'xlim');
count = length(y);
xs = repmat(xlims, count, 1);
ys = repmat(y(:), 1, 2);
h = plot(xs', ys', linetype{:});

end

author: ayonga

repo: frost-dev

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


text:
function createSureFitRawData(volPixSize)
% createSureFitRawData([volPixSize])
%
% volPixSize is in mm/pixel (defaults to [240/256 240/256 1.2])
%
% Converts Ifile data into the row-major ordered (8-bit uchar) data
% format requried by SureFit. This script is essentiall the createVAnatomy script
% with the output format modified slightly (in saveSureFitAnat()  -- which was
% saveVAnat())
%
% 00.01.20 RFD
%

if ~exist('volPixSize','var')
   volPixSize = [240/256 240/256 1.2];
   disp(['volPixSize defaulting to [ ' num2str(volPixSize,'%.4f ') ...
         	'].  I hope this is correct!']);
end

volume_pix_size = 1./volPixSize;

[fname, path] = uigetfile('*.*', 'Select one of the I-files...');

% load ifiles
%
disp('Loading I-files...');
img = makeCubeIfiles([path 'I'], 256, [1:124]);

disp('Finding optimal clip values...');
figure(99);
hist(img(:),100);
lowerClip = 0;
upperClip = 500;
answer = inputdlg({'lower clip value: ','upper clip value: '}, ...
   'Set intensity clip values', 1, ...
   {num2str(lowerClip),num2str(upperClip)}, 'on');
if ~isempty(answer)
	lowerClip = str2num(answer{1});
	upperClip = str2num(answer{2});
end
disp(['intensity clip values are: ' num2str(lowerClip) ', ' num2str(upperClip)]);

% Scale image values to be 0-255
% (if intensityClip <= 1, my scaleImage routine will clip the highest 'intensityClip'
% proportion of values, or, if intensityClip > 1, it clips it to the image value
% specified by intensityClip.
%
% I haven't been fully satisfied by just automatically clipping off the top,
% say 2% of intensities (this is what mrInitRet currently does to the inplanes).
% So, lately I've been looking at the histogram and picking intensityClip
% by hand, then viewing the images to see how well it worked.
disp('Clipping intensities...');
img(img<lowerClip) = lowerClip;
img(img>upperClip) = upperClip;
img = img-lowerClip;

% crop the image cube
% (you can skip this if you want, and just do the crop in mrGray)
figure(1);image(squeeze(img2(round(end/2),:,:))./clip.*255);colormap(gray(256));axis image;
figure(2);image(squeeze(img2(:,round(end/2),:))./clip.*255);colormap(gray(256));axis image;
figure(3);image(squeeze(img2(:,:,round(end/2)))./clip.*255);colormap(gray(256));axis image;
img2 = img2(20:200,30:240,1:124);
figure(1);image(squeeze(img(round(end/2),:,:))./clip.*255);colormap(gray(256));axis image;
figure(2);image(squeeze(img(:,round(end/2),:))./clip.*255);colormap(gray(256));axis image;
figure(3);image(squeeze(img(:,:,round(end/2)))./clip.*255);colormap(gray(256));axis image;

disp('Scaling to 0-255...');
img = round(img./upperClip*255);
disp('Original orientation:');
figure(99);
subplot(2,2,1);
image(squeeze(img(round(end/2),:,:)));colormap(gray(256));axis image;axis off;
subplot(2,2,2);
image(squeeze(img(:,round(end/2),:)));colormap(gray(256));axis image;axis off;
subplot(2,2,4);
image(squeeze(img(:,:,round(end/2))));colormap(gray(256));axis image;axis off;

% Reslice the data so that all voxels are cubic.
% The image currently has a field of view (fov) of 240mm in 256 voxels on the x and y axes
% and a fov of 148.800mm in 124 voxels on the z axis.

reslicedImg = resliceVoxels(img, 240, 240, 148.800);

% Reorient the data into the orientation expected by SureFit
% Basically, the data needs to be rotated 90 degrees around the anterior-posterior axis of the brain
% The image data is no 159x256x256

disp('Reorienting data...');

rotImg= zeros(planes,rows,cols);
for index=1:rows
   rotImg(:,index,:) = rot90(squeeze(reslicedImg(:,index,:)),3);
end

disp(['Volume dimensions are now:' int2str(size(rotImg))]);

disp('The reoriented data:');
figure(100);
subplot(2,2,1);
image(squeeze(rotImg(round(end/2),:,:)));colormap(gray(256));axis image;axis off;
subplot(2,2,2);
image(squeeze(rotImg(:,round(end/2),:)));colormap(gray(256));axis image;axis off;
subplot(2,2,4);
image(squeeze(rotImg(:,:,round(end/2))));colormap(gray(256));axis image;axis off;

% Save SureFitAnatomy
%
disp('Saving data in SureFit raw Format...');
path = saveSureFitAnat(rotImg);
save([path 'UnfoldParams'], 'volume_pix_size', 'lowerClip', 'upperClip');
disp(['SureFit raw data format saved to ' path]);
disp(['You should now run "python SUREFit/bin/Raw2Minc "' path size(rotImg)]);

return;

function path = saveSureFitAnat(imgCube, fileName)
% path = saveSureFitAnat(imgCube, [fileName])
% imgCube is a [rows x cols x planes] image array
% It must already be scaled to 0-255!
%
% fileName specifies the output file location (full path!)
% If omitted, then a save-file dialog box appears.
%
% The path to where the anatomy ends up is returned, in
% case you care.
%
%
% See Also: writeVolume, which does the same thing, but takes
% the input data in a different format.
%
% RFD

if ~exist('fileName', 'var')
   fileName = '';
end

% open file for writing (little-endian mode)
if isempty(fileName)
	[fname, path] = uiputfile('sureFitRawData.dat', 'Save SureFit raw data file...');
   fileName = [path fname];
else
   path = fileName;
end

vFile = fopen(fileName,'w','l');
if vFile<1
   while vFile<1;
   disp('Couldn''t open that file- please try saving it somewhere else.');
	[fname, path] = uiputfile('sureFitRawData.dat', 'Save SureFit raw data file...');
   fileName = [path fname];
   vFile = fopen(fileName,'w','l');
	end
end

% SureFit's Raw2Minc doesn't want any header, so nothing is written before the image data.

% Write data
count = fwrite(vFile, imgCube, 'uchar');
fclose(vFile);

return;


text:
function [r, v] = orb2eci(mu, oev)

% convert classical orbital elements to eci state vector

% input

%  mu     = gravitational constant (km**3/sec**2)
%  oev(1) = semimajor axis (kilometers)
%  oev(2) = orbital eccentricity (non-dimensional)
%           (0 <= eccentricity < 1)
%  oev(3) = orbital inclination (radians)
%           (0 <= inclination <= pi)
%  oev(4) = argument of perigee (radians)
%           (0 <= argument of perigee <= 2 pi)
%  oev(5) = right ascension of ascending node (radians)
%           (0 <= raan <= 2 pi)
%  oev(6) = true anomaly (radians)
%           (0 <= true anomaly <= 2 pi)

% output

%  r = eci position vector (kilometers)
%  v = eci velocity vector (kilometers/second)

% Orbital Mechanics with MATLAB

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

r = zeros(3, 1);
v = zeros(3, 1);

% unload orbital elements array

sma = oev(1);
ecc = oev(2);
inc = oev(3);
argper = oev(4);
raan = oev(5);
tanom = oev(6);

slr = sma * (1 - ecc * ecc);

rm = slr / (1 + ecc * cos(tanom));

arglat = argper + tanom;

sarglat = sin(arglat);
carglat = cos(arglat);

c4 = sqrt(mu / slr);
c5 = ecc * cos(argper) + carglat;
c6 = ecc * sin(argper) + sarglat;

sinc = sin(inc);
cinc = cos(inc);

sraan = sin(raan);
craan = cos(raan);

% position vector

r(1) = rm * (craan * carglat - sraan * cinc * sarglat);
r(2) = rm * (sraan * carglat + cinc * sarglat * craan);
r(3) = rm * sinc * sarglat;

% velocity vector

v(1) = -c4 * (craan * c6 + sraan * cinc * c5);
v(2) = -c4 * (sraan * c6 - craan * cinc * c5);
v(3) = c4 * c5 * sinc;


text:
function [nodes, error, C] = removeUnconnected(nodes,error,C)
% This function removes unconnected (dead) nodes from the DCS structure.

NumOfNodes = size(nodes,2);

i = 1;
while i<=NumOfNodes
          if (any(C(i,:))==0)

              C = [C(1:i-1,:); C(i+1:NumOfNodes,:);];
              C = [C(:,1:i-1)  C(:,i+1:NumOfNodes);];

              nodes = [nodes(:,1:i-1) nodes(:,i+1:NumOfNodes);];
                 error = [error(1,1:i-1) error(1,i+1:NumOfNodes);];

               NumOfNodes = NumOfNodes - 1;
               i = i-1;
          end
           i = i+1;
 end


text:
xk=[0,0]';
dk=-grad_obj2(xk);
alpha=linspace(0,0.1,50);
phi_alpha=zeros(1,50);
for i=1:50
    phi_alpha(i)=phi(alpha(i),xk,dk);
end

plot(alpha,phi_alpha);

repo: Optimization-Theory-and-Algorithm


text:
function [H,S,D]=sobi(X,n,p),
% program  by  A. Belouchrani and A. Cichocki
%
% Second
% Order
% Blind
% Identification
% SOBI
%**************************************************
% blind identification by joint diagonalization   *
% of correlation  matrices.                       *
%                       			  *
% ------------------------------------------------*
% THIS CODE ASSUMES TEMPORALLY CORRELATED SIGNALS *
% in estimating the cumulants                     *
% ------------------------------------------------*
%
% [H,S]=SOBI(X,m,p) produces a matrix H of dimension [m by n] and a matrix S
% of dimension [n by N] wich are respectively an estimate of the mixture matrix and
% an estimate of the source signals of the linear model from where the
% observation matrix X of dimension [m by N] derive.
% Note: > m: sensor number.
%       > n: source number by default m=n.
%       > N: sample number.
%	> p: number of correlation matrices to be diagonalized by default p=4.
%
% REFERENCES:
% A. Belouchrani, K. Abed-Meraim, J.-F. Cardoso, and E. Moulines, ``Second-order
%  blind separation of temporally correlated sources,'' in Proc. Int. Conf. on
%  Digital Sig. Proc., (Cyprus), pp. 346--351, 1993.
%
% A. Belouchrani and K. Abed-Meraim, ``Separation aveugle au second ordre de
%  sources correlees,'' in  Proc. Gretsi, (Juan-les-pins),
%  pp. 309--312, 1993.
%
%  A. Belouchrani, and A. Cichocki,
%  Robust whitening procedure in blind source separation context,
%  Electronics Letters, Vol. 36, No. 24, 2000, pp. 2050-2053.
%
%  A. Cichocki and S. Amari,
%  Adaptive Blind Signal and Image Processing, Wiley,  2002.

[m,N]=size(X);
if nargin==1,
 n=m; % source detection (hum...)

 p=4; % number of correlation matrices to be diagonalized
 end;
if nargin==2,
 p=4 ; % number of correlation matrices to be diagonalized
end;
pm=p*m; % for convenience
X=X-kron(mean(X')',ones(1,N)); % zero mean
%%%% whitening
% Rx=(X*X')/T;
% if m<n, %assumes white noise
%   [U,D]=eig(Rx); [puiss,k]=sort(diag(D));
%   ibl= sqrt(puiss(n-m+1:n)-mean(puiss(1:n-m)));
%    bl = ones(m,1) ./ ibl ;
%   BL=diag(bl)*U(1:n,k(n-m+1:n))';
%   IBL=U(1:n,k(n-m+1:n))*diag(ibl);
% else    %assumes no noise
%    IBL=sqrtm(Rx);
%    Q=inv(IBL);
% end;
% X=Q*X;
%prewitening based directly on SVD
[UU,S,VV]=svd(X',0);
Q= pinv(S)*VV';
X=Q*X;
%
%%%correlation matrices estimation
k=1;for u=1:m:pm, k=k+1; Rxp=X(:,k:N)*X(:,1:N-k+1)'/(N-k+1);
    M(:,u:u+m-1)=norm(Rxp,'fro')*Rxp; end;
%%%joint diagonalization
epsil=1/sqrt(N)/100; encore=1; V=eye(m);
while encore, encore=0;
 for p=1:m-1,
  for q=p+1:m,
   %%% Givens rotations
   g=[   M(p,p:m:pm)-M(q,q:m:pm)  ;
         M(p,q:m:pm)+M(q,p:m:pm)  ;
      i*(M(q,p:m:pm)-M(p,q:m:pm)) ];
	  [vcp,D] = eig(real(g*g')); [la,K]=sort(diag(D));
   angles=vcp(:,K(3));angles=sign(angles(1))*angles;
   c=sqrt(0.5+angles(1)/2);
   sr=0.5*(angles(2)-j*angles(3))/c; sc=conj(sr);
   oui = abs(sr)>epsil ;
   encore=encore | oui ;
   if oui , %%%update of the M and V matrices
    colp=M(:,p:m:pm);colq=M(:,q:m:pm);
    M(:,p:m:pm)=c*colp+sr*colq;M(:,q:m:pm)=c*colq-sc*colp;
    rowp=M(p,:);rowq=M(q,:);
    M(p,:)=c*rowp+sc*rowq;M(q,:)=c*rowq-sr*rowp;
    temp=V(:,p);
    V(:,p)=c*V(:,p)+sr*V(:,q);V(:,q)=c*V(:,q)-sc*temp;
   end%% if
  end%% q loop
 end%% p loop
end%% while
%%%estimation of the mixing matrix and source signals
Hf=pinv(Q)*V; %estimated mixing matrix
Sf=V'*X; %estimated sources

H = Hf(:,1:n);
S = Sf(1:n,:);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [p T hFig res] = rmCompareModelsGUI_paramTTest(M, plotFlag);
% Compare the distributions of pRF parameters between models in a compare
% models GUI.
%
%   [p T hFig df res] = rmCompareModelsGUI_paramTTest([M, plotFlag=1]);
%
% INPUTS:
%	M: rmCompareModelsGUI structure. See rmCompareModelsGUI_getData.
%	[Default: get from current figure]
%
%	plotFlag: flag to indicate whether to plot the results or just return
%	them. [default: 1, plot 'em]
%
% OUTPUTS:
%	p: [nModels x nModels x 7] matrix of p-values, for one-tailed T tests
%	between the parameters for each model. The 7 slices reflect the 7
%	parameters analyzed:
%		1) x0
%		2) y0
%		3) sigma (major)
%		4) polar angle
%		5) eccentricity
%		6) variance explained
%		7) beta coefficient for main pRF term.
%	So, for instance, p(2,1,3) reflects a comparison of the third parameter
%	(pRF size or sigma) for the test where model 2 is greater than model 1.
%
%	T: [nModels x nModels x 7] matrix of T-values, corresponding to the p-
%	values given in p.
%
%	hFig: if plotFlag==1, returns a handle to the plotted data. Otherwise,
%	returns empty.
%
%	res: further results struct, with the following fields:
%		df: degrees of freedom for the T tests
%
%		alpha: alpha threshold for the T tests, based on p=0.01 with
%		Bonferroni correction for multiple comparisons within these tests
%
%		H: matrix of 'results' indicating whether to accept (=0) or reject (=1)
%		the null hypothesis that the parameters come from the same
%		distribution. Same format as p and T.
%
%		CI_lo, CI_hi: 100 * (1-alpha)% confidence intervals for the true mean for
%		each comparison (same format as H). CI_lo is the lower bound, CI_hi
%		is the upper bound.
%
%		sd: pooled estimate of the population standard deviation for each
%		comparison (same format as H).
%
%
% ras, 05/2009.
if notDefined('M'),		M = get(gcf, 'UserData');			end
if notDefined('plotFlag'),	plotFlag = 1;					end
if ishandle(M),			M = get(M, 'UserData');				end

p = [];
T = [];
hFig = [];

fields = {'x0' 'y0' 'sigma' 'pol' 'ecc' 'varexp' 'beta'};

% let's pick an alpha value that takes into account multiple comparisons
% (not that it really matters -- in the initial iteration, we don't return
% these, and let the user decide their own thresholds based on the
% p-values):
alpha = 0.01 / (7 * M.nModels);
res.alpha = alpha;

%% loop across parameters
for n = 1:7
	f = fields{n};

	% get the param values for all models
	if n < 7
		data = reshape( [M.(f){:}], [M.nVoxels M.nModels] );
	else
		% n==7: beta values: need the specific beta for the main
		% effect
		for m = 1:M.nModels
			tmp{m} = [ M.beta{m}(:,1) ]';
		end
		data = reshape( [tmp{:}], [M.nVoxels M.nModels] );
	end

	%% loop across model pairs
	for ii = 1:M.nModels
		for jj = 1:M.nModels
			% run the one-sided T test
			[H P CI STATS] = ttest2(data(:,ii), data(:,jj), alpha, 'right');

			% record the main stat values
			p(ii,jj,n) = P;
			T(ii,jj,n) = STATS.tstat;

			% record other results
			res.df(ii,jj,n) = STATS.df;
			res.H(ii,jj,n) = H;
			res.CI_lo(ii,jj,n) = CI(1);
			res.CI_hi(ii,jj,n) = CI(2);
			res.sd(ii,jj,n) = STATS.sd;
		end
	end
end

%% plot if requested
if plotFlag==1
	nm = sprintf('Cross-Model T Tests %s', M.roi.name);
	h = figure('Color', 'w', 'Name', nm);

	for z = 1:7
		subplot(3, 3, z);

		% get the param values for all models
		f = fields{z};
		if z < 7
			data = reshape( [M.(f){:}], [M.nVoxels M.nModels] );
		else
			% ii==7: beta values: need the specific beta for the main
			% effect
			for m = 1:M.nModels
				tmp{m} = [ M.beta{m}(:,1) ]';
			end
			data = reshape( [tmp{:}], [M.nVoxels M.nModels] );
		end

		% compute the mean, SEM of the data
		Y = nanmean(data);
		E = nanstd(data); %  ./ sqrt(M.nVoxels - 1);

		starbar(Y, E, any(res.H(:,:,z)));
		set(gca, 'Box', 'off');  tuftify;
		xlabel('Model #', 'FontSize', 12);
		ylabel(f, 'FontSize', 12);

		title(f, 'FontSize', 14);
	end

	nm = sprintf('Cross-Model T Values %s', M.roi.name);
	h = figure('Color', 'w', 'Name', nm, 'Units', 'norm', ...
			  'Position', [.4 .3 .4 .4]);

	for z = 1:7
		subplot(3, 3, z);

		f = fields{z};

		drawXCorrMatrix( T(:,:,z), mrvMinmax(T) + [-.2 .2], 1 );

		% label each entry in the lower-left-hand plot
% 		if z==7
			axis on
			set(gca, 'Box', 'off', 'XTick', 1:M.nModels, ...
				'YTick', 1:1:M.nModels, 'FontSize', 9)
			xlabel('Model #', 'FontSize', 12);
			ylabel('Model #', 'FontSize', 12);
% 		end

		title( ['T Values: ' f], 'FontSize', 14 );
	end

end

return

repo: vistasoft


text:
function K = spikernel( k, d1, d2, ind1, ind2, kerparam)

% K = spikernel( k, d1, d2, ind1, ind2, kerparam)
% computes the spikernel between two sequences of spiking activity
%
% kerparam{1} = N; max subsequence lengths (the kernel compares subsequences of lengths 1 to N and returns a sum of kernels)
% kerparam{2} = lam; \lambda parameter from article
% kerparam{3} = mu; \mu parameter from article
% kerparam{4} = p; the q parameter in the article - the parameter that is used in the sum of kernels to weigh them differently
% kerparam{5} = bins; the number of bins used;          we need this to reshape the data correctly
% kerparam{6} = nneurons; the number of neurons;    we need this to reshape the data correctly
%
%   This code uses the spikernel function, defined in:
%   Spikernels: Embedding Spiking Neurons in Inner-Product Spaces
%   Lavi Shpigelman, Yoram Singer, Rony Paz and Eilon Vaadia
%   Advances in Neural Information Processing Systems (NIPS) 15
%   MIT Press, Cambridge, MA, 2003.
%
% for further details see Fspikernel.m

bins = kerparam{5};
nneurons = kerparam{6};

xx = get_x( d1, ind1);
yy = get_x( d2, ind2);
K = zeros( size( yy,1), size( xx,1));

for i = 1:size( xx, 1)
    for j = 1:size( yy,1)
        K( j, i) = Fspikernel( reshape( yy( j, :), bins, nneurons), reshape( xx( i,:), bins, nneurons), kerparam{1}, kerparam{2}, kerparam{3}, kerparam{4});
%        K( j, i) = K( i, j);
%       if bins >= 100 disp(sprintf('we''re at: (%d, %d)', i,j)); end   %<-- if it takes longer, print out, where we are
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore