text:
%%*********************************************************************
%% gdcomp: Compute gd = 1/td in Equation (15) of the paper:
%%
%% R.M. Freund, F. Ordonez, and K.C. Toh,
%% Behavioral measures and their correlation with IPM iteration counts
%% on semi-definite programming problems,
%% Mathematical Programming, 109 (2007), pp. 445--475.
%%
%% [gd,info,yfeas,Zfeas,blk2,At2,C2,b2] = gdcomp(blk,At,C,b,OPTIONS,solveyes);
%%
%% yfeas,Zfeas: a dual feasible pair when gd is finite.
%%              That is, if
%%              Aty = Atyfun(blk,At,[],[],yfeas);
%%              Rd = ops(C,'-',ops(Zfeas,'+',Aty));
%%              then
%%              ops(Rd,'norm') should be small.
%%
%%*********************************************************************

  function [gd,info,yfeas,Zfeas,blk2,At2,C2,b2] = gdcomp(blk,At,C,b,OPTIONS,solveyes);

  if (nargin < 6); solveyes = 1; end
  if (nargin < 5)
     OPTIONS = sqlparameters;
     OPTIONS.vers   = 1;
     OPTIONS.gaptol = 1e-10;
     OPTIONS.printlevel = 3;
  end
  if isempty(OPTIONS); OPTIONS = sqlparameters; end
  if ~isfield(OPTIONS,'solver'); OPTIONS.solver = 'HSDsqlp'; end
  if ~isfield(OPTIONS,'printlevel'); OPTIONS.printlevel = 3; end
  if ~iscell(C); tmp = C; clear C; C{1} = tmp; end
%%
  m = length(b);
  blk2 = blk;
  At2 = cell(size(blk,1),1);
  C2 = cell(size(blk,1),1);
  EE = cell(size(blk,1),1);
%%
%%
%%
  dd = zeros(1,m);
  alp = 0;
  beta = 0;
  for p = 1:size(blk,1)
     pblk = blk(p,:);
     n = sum(pblk{2});
     if strcmp(pblk{1},'s')
        C2{p,1} = sparse(n,n);
     else
        C2{p,1} = zeros(n,1);
     end
     dd = dd + sqrt(sum(At{p}.*At{p}));
     beta = beta + norm(C{p},'fro');
     alp = alp + sqrt(n);
  end
  alp  = 1./max(1,alp);
  beta = 1./max(1,beta);
  dd   = 1./max(1,dd);
%%
%% New multipliers in dual problem:
%% [v; tt; theta].
%%
   D = spdiags(dd',0,m,m);

   ss = 0; cc = 0; aa = zeros(1,m);
   exist_ublk = 0;
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      n = sum(pblk{2});
      if strcmp(pblk{1},'s')
         At2{p} = [At{p}*D, svec(pblk,alp*speye(n,n),1), -svec(pblk,beta*C{p},1)];
         ss = ss + n;
         cc = cc + trace(C{p});
         aa = aa + svec(pblk,speye(n),1)'*At{p};
         EE{p} = speye(n,n);
      elseif strcmp(pblk{1},'q')
         eq = zeros(n,1);
         idx1 = 1+[0,cumsum(pblk{2})];
         idx1 = idx1(1:length(idx1)-1);
         eq(idx1) = ones(length(idx1),1);
         At2{p} = [At{p}*D, 2*sparse(alp*eq), -sparse(beta*C{p})];
         ss = ss + 2*length(pblk{2});
         cc = cc + sum(C{p}(idx1));
         aa = aa + eq'*At{p};
         EE{p} = eq;
      elseif strcmp(pblk{1},'l')
         el = ones(n,1);
         At2{p} = [At{p}*D, sparse(alp*el), -sparse(beta*C{p})];
         ss = ss + n;
         cc = cc + el'*C{p};
         aa = aa + el'*At{p};
         EE{p} = el;
      elseif strcmp(pblk{1},'u')
         At2{p} = [At{p}*D, sparse(n,1), -sparse(beta*C{p})];
         exist_ublk = 1;
         EE{p} = sparse(n,1);
      end
   end
   aa = aa.*dd;
   cc = cc*beta;
%%
%% 4 additional inequality constraints in dual problem.
%%
   numblk = size(blk,1);
   blk2{numblk+1,1} = 'l'; blk2{numblk+1,2} = 4;
   C2{numblk+1,1}  = [1; 1; 0; 0];
   At2{numblk+1,1} = [-aa,         0,   cc;
		     zeros(1,m),   0,   beta;
		     zeros(1,m),  alp, -beta
                     zeros(1,m), -alp,  0];
   At2{numblk+1} = sparse(At2{numblk+1});
   b2 = [zeros(m,1); alp; 0];
%%
%% Solve SDP
%%
   gd = []; info = []; yfeas = []; Zfeas = [];
   if (solveyes)
      if strcmp(OPTIONS.solver,'sqlp')
         [X0,y0,Z0] = infeaspt(blk2,At2,C2,b2,2,100);
         [obj,X,y,Z,info] = sqlp(blk2,At2,C2,b2,OPTIONS,X0,y0,Z0);
      else
         [obj,X,y,Z,info] = HSDsqlp(blk2,At2,C2,b2,OPTIONS);
      end
      tt = alp*abs(y(m+1)); theta = beta*abs(y(m+2));
      yfeas = D*y(1:m)/theta;
      Zfeas = ops(ops(Z(1:numblk),'+',EE,tt),'/',theta);
      %%
      if (obj(2) > 0) | (abs(obj(2)) < 1e-8)
         gd = 1/abs(obj(2));
      elseif (obj(1) > 0)
         gd = 1/obj(1);
      else
         gd = 1/exp(mean(log(abs(obj))));
      end
      err = max(info.dimacs([1,3,6]));
      if (OPTIONS.printlevel)
         fprintf('\n ******** gd = %3.2e, err = %3.1e\n',gd,err);
         if (err > 1e-6);
            fprintf('\n----------------------------------------------------')
            fprintf('\n gd problem is not solved to sufficient accuracy');
            fprintf('\n----------------------------------------------------\n')
         end
      end
   end
%%*********************************************************************

repo: REST


text:
% usage demo for

% Masked k-order statistic filter for double data
% -------------------------------------------------
% by Fabio Bellavia (fbellavia@unipa.it),
% refer to: F. Bellavia, D. Tegolo, C. Valenti,
% "Improving Harris corner selection strategy",
% IET Computer Vision 5(2), 2011.
% Only for academic or other non-commercial purposes.

% see kordstatfilt2.m for more details

% compile the mex file
% tested with Matlab R2012a x32 on Ubuntu 12.04
mex ordstatfilt2.c

% input data
im=double(rgb2gray(imread('donkey.jpg')));
figure;
imshow(im,[]);

% circular kernel
ker=fspecial('disk',15)>0;

% median filter
idx=sum(ker(:))/2+0.5;
res=kordstatfilt2(im,ker,idx);
figure;
imshow(res,[]);

author: Sable

repo: mcbench-benchmarks


text:
function [ n_data, x, fx ] = goodwin_values ( n_data )

%*****************************************************************************80
%
%% GOODWIN_VALUES returns some values of the Goodwin and Staton function.
%
%  Discussion:
%
%    The function is defined by:
%
%      GOODWIN(x) = Integral ( 0 <= t < infinity ) exp ( -t^2 ) / ( t + x ) dt
%
%    The data was reported by McLeod.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%      special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  fx_vec = [ ...
      0.59531540040441651584E+01, ...
      0.45769601268624494109E+01, ...
      0.32288921331902217638E+01, ...
      0.19746110873568719362E+01, ...
      0.96356046208697728563E+00, ...
      0.60513365250334458174E+00, ...
      0.51305506459532198016E+00, ...
      0.44598602820946133091E+00, ...
      0.37344458206879749357E+00, ...
      0.35433592884953063055E+00, ...
      0.33712156518881920994E+00, ...
      0.29436170729362979176E+00, ...
      0.25193499644897222840E+00, ...
      0.22028778222123939276E+00, ...
      0.19575258237698917033E+00, ...
      0.17616303166670699424E+00, ...
      0.16015469479664778673E+00, ...
      0.14096116876193391066E+00, ...
      0.13554987191049066274E+00, ...
      0.11751605060085098084E+00 ];

  x_vec = [ ...
      0.0019531250E+00, ...
      0.0078125000E+00, ...
      0.0312500000E+00, ...
      0.1250000000E+00, ...
      0.5000000000E+00, ...
      1.0000000000E+00, ...
      1.2500000000E+00, ...
      1.5000000000E+00, ...
      1.8750000000E+00, ...
      2.0000000000E+00, ...
      2.1250000000E+00, ...
      2.5000000000E+00, ...
      3.0000000000E+00, ...
      3.5000000000E+00, ...
      4.0000000000E+00, ...
      4.5000000000E+00, ...
      5.0000000000E+00, ...
      5.7500000000E+00, ...
      6.0000000000E+00, ...
      7.0000000000E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function model = loss_pyramid(h_loss_func, pyra, model, fg_box, ...
                              bg_boxes, min_fg_overlap, max_bg_overlap)
% Computes a pyramid of loss function values for each top-level
% rule in the grammar.
%   model = loss_pyramid(h_loss_func, pyra, model, fg_box, ...
%                        bg_boxes, min_fg_overlap, max_bg_overlap)
%
%   These loss values are used for computing the loss adjusted inference:
%     \max_{s \in S(x)} w \dot \psi(x,s) + L_margin(y,s)
%   The set of valid outputs S(x) is enforced by making L(y,s) = -inf for
%   some values of s, which prevents them from being selected in the
%   maximization.
%
% Return value
%   model           Model augmented to store the computed loss pyramids
%
% Arguments
%   h_loss_func     Handle to loss function
%   model           Model
%                   (augmented with DP tables from gdetect_dp.m)
%   pyra            Feature pyramid
%                   (augmented with overlaps from gdetect_pos_prepare.m)
%   fg_box          Selected foreground bounding box index
%   bg_boxes        Indices of non-selected bounding boxes in image
%   min_fg_overlap  Minimum required amount of overlap with fg box
%   max_bg_overlap  Maximum allowed amount of overlap with bg bounding boxes

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

num_bg_boxes = length(bg_boxes);

% For each model component
for comp = 1:length(model.rules{model.start})
  % For each feature pyramid level
  for level = 1:pyra.num_levels
    if pyra.valid_levels(level)
      % Assign loss for root locations based on the selected foreground box
      o = pyra.overlaps(comp).box(fg_box).o{level};
      losses = h_loss_func(o);
      model.rules{model.start}(comp).loss{level} = losses;

      % Require at least some overlap with the foreground bounding box
      % Rationale:
      %  In an image with multiple objects, this constraint encourages a
      %  diverse set of false positives (otherwise, they will tend to come
      %  from the same high-scoring / low-overlapping region of the image
      %  -- i.e. somewhere in the background)
      I = find(o < min_fg_overlap);
      model.rules{model.start}(comp).loss{level}(I) = -inf;

      % Mark root locations that have too much overlap with background boxes
      % as invalid
      % Rationale:
      %   We don't want to select detections of other foreground objects
      %   in the image as false positives (i.e., no true positive should
      %   be allowed to be used as a false positive)
      for b = 1:num_bg_boxes
        o = pyra.overlaps(comp).box(bg_boxes(b)).o{level};
        inds = find(o >= max_bg_overlap);
        model.rules{model.start}(comp).loss{level}(inds) = -inf;
      end
    else
      model.rules{model.start}(comp).loss{level} = 0;
    end
  end
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function test_suite=test_moxunit_util_regexp_matches()
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function c=randchar(length)
    c=char(ceil(rand(1,length)*26+64));

function test_moxunit_util_regexp_matches_basics
    for needle_length=1:3

        %  make string to match
        needle=randchar(needle_length);
        assert_regexp_matches(needle,needle)

        for pat_length=1:4
            % make pattern that matches or does not match needle
            pat=randchar(pat_length);
            assert_regexp_matches([needle pat],needle)
            assert_regexp_matches([pat needle],needle)

            if needle_length==pat_length
                matcher=@assert_regexp_matches;
            else
                matcher=@assert_regexp_not_matches;
            end
            matcher(needle,sprintf('^%s$',...
                                        repmat('.',1,pat_length)))

            while true
                different_haystack=randchar(pat_length);
                match=bsxfun(@eq,different_haystack,needle');
                if all(~match(:))
                    break;
                end
            end

            assert_regexp_not_matches(needle,different_haystack)
        end
    end

function test_moxunit_util_regexp_matches_exceptions()
    aet=@(varargin)assertExceptionThrown(@()...
                        moxunit_util_regexp_matches(varargin{:}),'');
    bad_inputs={{[],''},... % first argument not a string
                 {'',[]},... % second argument not a string
                 };
    for k=1:numel(bad_inputs)
        args=bad_inputs{k};
        aet(args{:});
    end

function assert_regexp_matches(needle, pat)
    assertTrue(moxunit_util_regexp_matches(needle,pat));

function assert_regexp_not_matches(needle, pat)
    assertFalse(moxunit_util_regexp_matches(needle,pat));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
data_filename = ['generated_scenes.mat'];

%% params to be set
data_folder = ['..', filesep, '0-data'];
input_folder = [data_folder, filesep, '4-pydata']
labelset_folder = [data_folder, filesep, 'labelsets'];
output_folder = [data_folder, filesep, '5-generated_scenes'];
output_obj_folder = [output_folder, filesep, 'recon'];
output_img_folder = [output_folder, filesep, 'images'];
objfolder_path = [data_folder, filesep, 'SUNCG', filesep, 'object'];
labelset_filename = [labelset_folder,filesep,'bedroomlabelset.mat'];
objperlabel_filename = [data_folder,filesep,'object_per_label',filesep,'bedroom_objset_per_label.mat'];

addpath(genpath('..\3-datapreparation'));

if ~exist(output_folder,'dir')
   mkdir(output_folder);
end
if ~exist(output_obj_folder, 'dir')
	mkdir(output_obj_folder);
end
if ~exist(output_img_folder, 'dir')
	mkdir(output_img_folder);
end

repo: GRAINS


text:
% TEST_WAVELET_FACTORY_2D Test case for WAVELET_FACTORY_2D
%
% See also
%   WAVELET_FACTORY_2D
classdef test_wavelet_factory_2d < matlab.unittest.TestCase
    methods(Test)

        function testWithNoOptions(testcase)
            %% define
            x = rand(64, 64);
            [Wop, filters] = wavelet_factory_2d(size(x));

            %% check number of high pass filter
            expected = 3;
            actual = numel(Wop);

            %% assert
            testcase.assertEqual(expected, actual);

        end

          function testWithRandomOptions(testcase)
            for i = 1:32
                %%
                sz = 1 + floor(128*rand(1,2));
                %% define
                white_list_filt = {'Q', 'J', 'L'};
                type={'d','d','d'};
                values={{1,2,4,8,16},{0,1},{1,10}};
                options_f = generate_random_options(white_list_filt,type,values);

                white_list_scat = { 'oversampling','M'};
                type={'d','i'};
                values={{0,1},{1,3}};

                options_s = generate_random_options(white_list_scat,type,values);

                x = rand(64, 64);
                [Wop, filters] = wavelet_factory_2d(size(x),options_f,options_s);

                actual = length(Wop);
                expected = options_s.M;
                %% assert
                testcase.assertEqual(expected, actual);
            end

        end

        function testWithRandomSize(testcase)
            for i = 1:32
                %%
                sz = 1 + floor(128*rand(1,2));
                %% define
                x = rand(sz);
                [Wop, filters] = wavelet_factory_2d(size(x));

                % Test
                expected = 3;
                actual = numel(Wop);

                %% assert
                testcase.assertEqual(expected, actual);

            end
            for i = 1:32
                %%
                sz = 1 + floor(2*rand(1,2));
                %% define
                x = rand(sz);
                [Wop, filters] = wavelet_factory_2d(size(x));

                % Test
                expected = 3;
                actual = numel(Wop);

                %% assert
                testcase.assertEqual(expected, actual);

            end
    end
    end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [outputStruct]=febioStructTemplate_v2p5(varargin)

switch nargin
    case 0
        inputStruct=[];
    case 1
        inputStruct=varargin{1};
    otherwise
        error('Wrong number of input arguments');
end

%% Set febio_spec version
febio_spec.ATTR.version='2.5';

%% Module section
febio_spec.Module.ATTR.type='solid'; %Use default set

%% Control section

febio_spec.Control.analysis.ATTR.type='static';
% febio_spec.Control.title='temp';
febio_spec.Control.time_steps=10;
febio_spec.Control.step_size=0.1;
febio_spec.Control.dtol=0.001;
febio_spec.Control.etol=0.01;
febio_spec.Control.rtol=0;
febio_spec.Control.lstol=0.9;
febio_spec.Control.time_stepper.dtmin=febio_spec.Control.step_size/3;
febio_spec.Control.time_stepper.dtmax=febio_spec.Control.step_size*3; %Constant
% febio_spec.Control.time_stepper.dtmax.ATTR.lc=1; %Load curve based
% febio_spec.Control.time_stepper.dtmax.VAL=0.1; %Load curve based
febio_spec.Control.time_stepper.max_retries=5;
febio_spec.Control.time_stepper.opt_iter=10;
febio_spec.Control.max_refs=15;
febio_spec.Control.max_ups=10;
febio_spec.Control.optimize_bw=0;
% febio_spec.Control.restart.ATTR.file='restartDumpFile.dmp'; %On with file specified
% febio_spec.Control.restart.VAL=0; %Off
febio_spec.Control.plot_level='PLOT_MAJOR_ITRS';
febio_spec.Control.plot_range=[0,-1];
febio_spec.Control.plot_stride=1;
febio_spec.Control.plot_zero_state=0;
febio_spec.Control.cmax=1e5;
febio_spec.Control.print_level='PRINT_MINOR_ITRS';
febio_spec.Control.min_residual=1e-20;
% febio_spec.Control.integration='N/A';
febio_spec.Control.output_level='OUTPUT_MAJOR_ITRS';

%% Globals section

febio_spec.Globals.Constants.T=0;
febio_spec.Globals.Constants.R=0;
febio_spec.Globals.Constants.Fc=0;

%% Material section
% febio_spec.Material.material{1}.ATTR.type='Ogden';
% febio_spec.Material.material{1}.ATTR.id=1;
% febio_spec.Material.material{1}.c1=1e-3;
% febio_spec.Material.material{1}.m1=6;
% febio_spec.Material.material{1}.c2=1e-3;
% febio_spec.Material.material{1}.m2=-6;
% febio_spec.Material.material{1}.k=1e-1;

%% LoadData section

febio_spec.LoadData.loadcurve{1}.ATTR.id=1;
febio_spec.LoadData.loadcurve{1}.ATTR.type='linear';
febio_spec.LoadData.loadcurve{1}.point.VAL=[0 0; 1 1];

%% Output section

% Plot file
febio_spec.Output.plotfile.ATTR.type='febio';
febio_spec.Output.plotfile.var{1}.ATTR.type='displacement';
febio_spec.Output.plotfile.var{2}.ATTR.type='stress';
febio_spec.Output.plotfile.var{3}.ATTR.type='strain energy density';
febio_spec.Output.plotfile.var{4}.ATTR.type='relative volume';
febio_spec.Output.plotfile.var{5}.ATTR.type='reaction forces';
febio_spec.Output.plotfile.var{6}.ATTR.type='contact pressure';

%% Fill in missing if input structure is provided

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function y = polyval( p, x )

%   Disciplined convex/geometric programming information for POLYVAL:
%      POLYVAL can be used with CVX variables in two ways:
%      -- If P is constant and X is a variable, then POLYVAL(P,X)
%         represents a polynomial function of X, p(X). If p(X) is affine
%         then the basic  DCP or DCP rules for sums and products apply.
%         If P is nonlinear then:
%         -- for DCPs, the polynomial must be convex or concave, and X
%            must be affine (such polynomials cannot be monotonic).
%         -- for DGPs, the elements of P must be nonnegative, in which
%            case p(X) is nondecreasing, so X must be log-convex.
%      -- If X is a constant and P is a variable, then POLYVAL(P,X)
%         performs a weighted sum of the elements of P. The weights are,
%         of course, [X^N,X^(N-1),...,X,1]. As in the affine case above,
%         then standard DCP or DGP rules for sums and products apply.

sp = size( p );
if isempty( p ),
    p = zeros( 1, 0 );
elseif length( sp ) > 2 || ~any( sp == 1 ),
    error( 'First argument must be a vector.' );
end
sx = size(x);

persistent remap
if isempty( remap ),
    remap_1 = cvx_remap( 'real-affine' );
    remap_2 = cvx_remap( 'log-convex' );
    remap = remap_1 + 2 * remap_2;
end

if cvx_isconstant( p ),
    p = cvx_constant( p );
    if cvx_isconstant( x ),
        y = cvx( polyval( p, cvx_constant( x ) ) );
        return
    end
    if any( isinf( p ) | isnan( p ) ),
        error( 'Inf and NaN not accepted here.' );
    end
    ndxs = find( p );
    if isempty( ndxs ),
        y = zeros( sx );
        return
    end
    for k = ndxs(:)',
        pt = p( k : end );
        if ~any( isinf( pt ./ pt(1) ) ),
            p = pt;
            break;
        end
    end
    n = length( p );
    switch length( p ),
        case 0,
            % Zero
            y = zeros(sx);
        case 1,
            % Constant
            y = cvx( p(1) * ones(sx) );
        case 2,
            % Affine
            y = p(1) * x + p(2);
        otherwise,
            vu = remap(cvx_classify(x));
            p = reshape( p, n, 1 );
            t0 = vu == 0;
            if any( t0 ),
                error( 'Disciplined convex programming error:\n    Illegal operation: polyval( p, {%s} ).', cvx_class( cvx_subsref( x, t0 ), false, true ) );
            end
            t1 = vu == 1;
            if any( t1 ),
                pd = roots( (n-2:-1:1).*(n-1:-1:2).*p(1:end-2,:)' );
                pd = pd( imag(pd) == 0 );
                if ~isempty( pd ),
                    pr = diff( [ 0, find(diff(sort(pd))~=0), length(pd) ] );
                    if any( rem( pr, 2 ) ),
                        error( 'Polynomials in DCPs must be affine, convex, or concave.' );
                    end
                end
            end
            t2 = vu == 2;
            if any( t2 ) && any( p < 0 ),
                error( 'Polynomials in DGPs must have nonnegative coefficients.' );
            end
            if any( t2 ),
                if all( t2 ),
                    y = x;
                else
                    y = cvx_subsref( x, t2 );
                end
                ny = numel(y);
                y = reshape( y, ny, 1 );
                q = find( p(1:end-2) );
                y = ( ( y * ones(1,length(q)) ) .^ ( ones(ny,1) * (n-q)' ) ) * p(q,:) + p(end-1) * y + p(end);
                if all( t2 ),
                    y = reshape( y, sx );
                else
                    y = cvx_subsasgn( x, t2, y );
                end
            end
            if all( t1 ),
                y = poly_env( p, x );
            elseif any( t1 ),
                y = cvx_subsasgn( y, t1, poly_env( p, cvx_subsref( x, t1 ) ) );
            end
    end

elseif cvx_isconstant( x ),

    n = length( p );
    [ ii, jj, vv ] = find( p );
    jj = ii + jj - 1;
    nv = length(vv);
    nx = prod( sx );
    y = reshape( x, nx, 1 ) * ones( 1, nv );
    y = y .^ ( ones( nx, 1 ) * (n-jj(:))' );
    y = reshape( y * reshape( vv, nv, 1 ), sx );

else

    error( 'At least one of the arguments must be constant.' );

end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function rts_updateConnStatusLbl(connected, hLbl, vesselName)
%rts_updateConnStatusLbl Summary of this function goes here
%   Detailed explanation goes here

    connStr = 'CONNECTED';
    disconnStr = 'DISCONNECTED';
    if(connected)
        preposition = ' to ';
        statusStr = formConnectedStatusString(connStr, vesselName, preposition);
        statusColor = 'green';
    else
        statusStr = formDisconnectedStatusString(disconnStr);
        statusColor = 'red';
    end

    if(ishandle(hLbl))
        set(hLbl, 'String', statusStr);
        set(hLbl, 'ForegroundColor', statusColor);
    end
end

function statusStr = formConnectedStatusString(status, vesselName, preposition)
%     statusStr = [status, preposition, vesselName];
    statusStr = status;
end

function statusStr = formDisconnectedStatusString(status)
    statusStr = status;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [gx] = g_rbf(x,P,u,in)

centers = in.centers + exp(P(1));
sig = exp(P(2));

N = length(in.centers);
Xrbf = zeros(length(in.grid),N);
for i=1:N
    Xrbf(:,i) = exp(-0.5*(centers(i)-in.grid).^2./sig);
    Xrbf(:,i) = Xrbf(:,i)./sum(Xrbf(:,i));
end

gx = zeros(N,1);
for i=1:N
    y = Xrbf(:,i);
    X = Xrbf(:,setdiff(1:N,i));
    P0 = eye(length(in.grid)) - X*pinv(X'*X)*X';
    err = y'*P0*y;
    gx(i) = sum(err);
end

gx = gx + P(3);

if ~in.corr
    gx = Xrbf;
end

author: MBB-team

repo: VBA-toolbox


text:
function str = verbatim
%VERBATIM  Get the text that appears in the next comment block.
%  Returns the text of the first comment block following the call.  The
%  block comment delimiters, %{ and %}, must appear on lines by themselves
%  (optionally preceded by white space).
%
%  If you want a final end-of-line character then leave a blank line before
%  the %}.
%
%  If both comment delimiters are preceded by the same white space (same
%  combination of spaces and tabs) then that white space will be deleted
%  (if possible) from the beginning of each line of the commented text.
%  This is so the whole block can be indented.
%
%  Example,
%
%      str = verbatim;
%          %{
%          This is the text
%          that will be returned by verbatim.
%          %}
%
%  VERBATIM can only be used in an m-file.

% Version: 1.1, 22 March 2009
% Author:  Douglas M. Schwarz
% Email:   dmschwarz=ieee*org, dmschwarz=urgrad*rochester*edu
% Real_email = regexprep(Email,{'=','*'},{'@','.'})

% Get the function call stack.
[dbs,thisWorkspace] = dbstack('-completenames');
assert(length(dbs) > 1,'VERBATIM must be called from an M-file.')
dbs = dbs(thisWorkspace + 1);
lines = repmat({''},1,100);

% Open the file.
fid = fopen(dbs.file);

try
	% Skip lines up to the current line in the calling function.
	textscan(fid,'%*s',0,'HeaderLines',dbs.line,'Delimiter','');

	% Read lines until one begins with '%{'.
	line = '';
	while ~isequal(line,-1) && isempty(regexp(line,'^\s*%{','once'))
		line = fgetl(fid);
	end

	% Read and save lines until one begins with '%}'.  Leave first cell
	% empty.
	k = 1;
	while true
		k = k + 1;
		lines{k} = fgetl(fid);
		if isequal(lines{k},-1) || ...
				~isempty(regexp(lines{k},'^\s*%}','once'))
			break
		end
	end

	% Close the file.
	fclose(fid);

catch err
	% Close the file and rethrow the error.
	fclose(fid);
	rethrow(err)
end

% If white space preceeding '%{' and '%}' is the same then delete it from
% the beginning of each line of the commented text so you can have indented
% code.
white_space1 = regexp(line,'^\s*','match','once');
white_space2 = regexp(lines{k},'^\s*','match','once');
if strcmp(white_space1,white_space2)
	lines = regexprep(lines,['^',white_space1],'');
end

% Construct the output string.
str = [sprintf('%s\n',lines{2:k-2}),lines{k-1}];


text:
function yearly_tmin_correction(filenameout)

% this program is design to adjust the daily Tmin generated by WG using
% new yearly Tmin series produced by FFT.

% load the yearly averaged Tmin after FFT
load('Pnew_tmin');
tmin_FFT=Pnew_tmin';

% load WG generated daily Tmin
load(filenameout);
[n,m]=size(gTmin);
gTmin=gTmin';
tmin_WG=reshape(gTmin,[],1);

n=length(tmin_FFT); % years
m=length(tmin_WG);  % days

% calculate yearly Tmin generated by WG, namely (Y)
j=1;
Z=zeros(n,1);
for i=1:365:m
    Z(j,1)=mean(tmin_WG(i:i+365-1,1));
    j=j+1;
end

% calculate the ratio of yearly Tmin between those after FFT and
% generated by weather generator
for i=1:n
    tmin_ratio(i,1)=tmin_FFT(i,1)-Z(i,1);
end

% extend the yearly Tmin ratio to daily scale,the data in each
% year are the same
tmin_extent=zeros(m,1);
j=1;
for i=1:365:m
    tmin_extent(i:i+365-1,1)=tmin_ratio(j,1);
    j=j+1;
end

% adjust the daily Tmin generated by WG using above ratios
tmin_adjust=zeros(size(tmin_WG));
for i=1:m
    tmin_adjust(i,1)=tmin_WG(i,1)+tmin_extent(i,1);
end
yearly_corrected_tmin=tmin_adjust;
save('yearly_corrected_tmin','yearly_corrected_tmin')


text:
function fg = fgRead(fgFile)
% Read fiber group file to create a FG structure.
%
%  fgRead([fgFile=mrvSelectFile])
%
% INPUTS:
%   fgFile - A fiber group on file. C
%          - Compatible file formats are:
%            a. VISTASOFT .mat and .pdb
%            b. MRTrix .tck
%            c. TrackVis .trk
%
% WEB RESOURCES:
%   mrvBrowseSVN('fgRead');
%
% EXAMPLE:
%   >> fg = fgRead('fgName.mat');
%   >> fg = fgRead('fgName.pdb')
%   >> fg = fgRead('fgName.tck')
%   >> fg = fgRead('feName.trk')
%
% See Also:
%   fgWrite.m
%
%
% (C) Stanford VISTA, 2017 | Pestilli Franco

%% Check inputs
% Check for name variable and use fg.name if empty
if ~exist('fgFile','var') || isempty(fgFile)
    fgFile = mrvSelectFile('r',{'pdb','mat','tck'},'Select FG File');
end

% Get the file type
[~,~,fileType] = fileparts(fgFile);

%% Read the file
switch fileType
    case '.pdb'
        fg = mtrImportFibers( fgFile );
    case '.mat'
        fg = dtiLoadFiberGroup( fgFile );
    case '.tck'
        fg = dtiImportFibersMrtrix( fgFile );
    case '.trk'
        fg = read_trk_to_fg( fgFile );
    otherwise
        error('[%s] Cannot parse file type for the tractogram.',mfilename)
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%kcmaxlik_printsig 'Print a Maximum Likelihood Signature '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros cmaxlik_printsig.pane file
%
% Parameters:
% InputFile: isig 'Input Signature', required: 'input signature file'
% MultiChoice: omode 'Output Mode', default: 1: 'output mode toggle'
%    Choices are:
%   1: 'Text'
%   2: 'LaTeX'
%   3: 'HTML'
%   4: 'CSV'
% OutputFile: o 'Output', optional: 'output signature info'
%
% Example: o = kcmaxlik_printsig(isig, {'isig','';'omode',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% cmaxlik_printsig - Print a Maximum Likelihood Signature
%
%  DESCRIPTION
% This operator prints out the contents of a Maximum Likelihood Signature file created with "cmaxlik_signature". It expects a valid Maximum Likelihood signature file (passed with the parameter [-isig]), so its input must be either the output of a \fIcmaxlik_signature\fP or the output of \fIkappend\fP which appended several signatures together in the time dimension. Its output (passed with the parameter [-o], or redirected to kstdout) can be formatted so it can be used in reports, etc.
% The output of this program can be in four formats, depending on the value of the parameter [-omode] :
% The Text (ASCII) format option ([-omode 1]) formats the output so it can be viewed with simple text editors or in kstdout.
% The LaTeX format option ([-omode 2]) formats the output so it can be included in a LaTeX file with the \\input command.
% The HTML format option ([-omode 3]) formats the output in a HTML table so it can be read into a HTML file.
% The CSV format option ([-omode 4]) formats the output in a Comma-delimited format so it can be read by spreadsheets, etc.
%
%
%
%  EXAMPLES
% All examples for the Classify toolbox are listed on the Classify Toolbox Manual. For an example of this routine, please see the example workspace Classify:workspaces:MAXLIK-Signatures.
%
%  "SEE ALSO"
% cmaxlik_signature
%
%  RESTRICTIONS
% Expect the input to be a valid signature file for the Maximum Likelihood classifier.
%
%  REFERENCES
% All references for the Classify toolbox are listed on the Classify Toolbox Manual.
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

function varargout = kcmaxlik_printsig(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kcmaxlik_printsig(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'isig', '__input';'omode', 1;'o', '__output'};
maxval={0,0,1};
minval={0,0,1};
istoggle=[0,0,1];
was_set=istoggle * 0;
paramtype={'InputFile','MultiChoice','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'cmaxlik_printsig"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%BagOfWords Bag of words class
%
% The BagOfWords class holds sets of features for a number of images and
% supports image retrieval by comparing new images with those in the 'bag'.
%
% Methods::
% isword        Return all features assigned to word
% occurrences   Return number of occurrences of word
% remove_stop   Remove stop words
% wordvector    Return word frequency vector
% wordfreq      Return words and their frequencies
% similarity    Compare two word bags
% contains      List the images that contain a word
% exemplars     Display examples of word support regions
% display       Display the parameters of the bag of words
% char          Convert the parameters of the bag of words to a string
%
% Properties::
% K         The number of clusters specified
% nstop     The number of stop words specified
% nimages   The number of images in the bag
%
% Reference::
%
% J.Sivic and A.Zisserman,
% "Video Google: a text retrieval approach to object matching in videos",
% in Proc. Ninth IEEE Int. Conf. on Computer Vision, pp.1470-1477, Oct. 2003.
%
% See also PointFeature.

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

classdef BagOfWords < handle
    properties
        features    % vector of all features in the bag (PointFeature class)

        K       % number of clusters
        C       % cluster centres (NDxNW)
        words   % vector of word indices (NW)

        nstop     % number of stop words
        stopwords  % list of stop words
        map        % maps word index with stop words to word index without stop words

        nimages       % number of images (NI)

        wv          % cached word vectors
    end

    methods
        function bag = BagOfWords(sf, a1)
        %BagOfWords.BagOfWords Create a BagOfWords object
        %
        % B = BagOfWords(F, K) is a new bag of words created from the feature
        % vector F and with K words.  F can also be a cell array, as produced
        % by ISURF() for an image sequence.
        %
        % The features are sorted into K clusters and each cluster is termed
        % a visual word.
        %
        % B = BagOfWords(F, B2) is a new bag of words created from the feature
        % vector F but clustered to the words (and stop words) from the existing
        % bag B2.
        %
        % Notes::
        % - Uses the MEX function vl_kmeans to perform clustering (vlfeat.org).
        %
        % See also PointFeature, ISURF.

            % save the feature vector
            if iscell(sf)
                bag.features = [sf{:}];
            else
                bag.features = sf;
            end
            bag.nimages = max([bag.features.image_id]);

            if isnumeric(a1)
                K = a1;
                % do the clustering
                [bag.C,L] = vl_kmeans([bag.features.descriptor], K, ...
                    'verbose', 'algorithm', 'elkan');

                bag.K = K;
                bag.words = double(L);

            elseif isa(a1, 'BagOfWords')
                oldbag = a1;

                % cluster using number of words from old bag
                bag.words = closest([bag.features.descriptor], oldbag.C);

                if oldbag.stopwords > 0
                    % remove stopwords as per the original bag file
                    bag.K = oldbag.K;
                    bag.stopwords = oldbag.stopwords;
                    k = find(ismember(bag.words, oldbag.stopwords));

                    fprintf('Removing %d features associated with stop words\n', length(k));

                    bag.words(k) = [];
                    bag.words = oldbag.map(bag.words);
                    bag.features(k) = [];
                end

                bag.compute_wv(oldbag);

            end
        end

        function f = isword(bag, words)
        %BagOfWords.isword Features from words
        %
        % F = B.isword(W) is a vector of feature objects that are assigned to any of
        % the word W.  If W is a vector of words the result is a vector of features
        % assigned to all the words in W.
            k = ismember(bag.words, words);
            f = bag.features(k);
        end

        function n = occurrence(bag, word)
        %BagOfWords.occurrence Word occurrence
        %
        % N = B.occurrence(W) is the number of occurrences of the word W across
        % all features in the bag.
            n = sum(bag.words == word);
        end

        function [all2, S] = remove_stop(bag, nstop)
        %BagOfWords.remove_stop Remove stop words
        %
        % B.remove_stop(N) removes the N most frequent words (the stop words)
        % from the bag.  All remaining words are renumbered so that the word
        % labels are consecutive.

            [w,f] = count_unique(bag.words);
            [f,i] = sort(f, 'descend');
            bag.stopwords = w(i(1:nstop));

            % remove all features that are stop words from L and all
            k = find(ismember(bag.words, bag.stopwords));

            fprintf('Removing %d features associated with %d most frequent words\n', ...
                length(k), nstop);

            % fix the labels
            b = zeros(1,length(bag.words));
            b(bag.stopwords) = 1;
            bag.map = [1:length(bag.words)] - cumsum(b);

            bag.words(k) = [];
            bag.words = bag.map(bag.words);
            bag.features(k) = [];

        end

        function wv = wordvector(bag, k)
        %BagOfWords.wordvector  Word frequency vector
        %
        % WF = B.wordvector(J) is the word frequency vector for the J'th image
        % in the bag.  The vector is Kx1 and the angle between any two WFVs is
        % an indication of image similarity.
        %
        % Notes::
        % - The word vector is expensive to compute so a lazy evaluation is
        %   performed on the first call to this function
            if isempty(bag.wv)
                bag.compute_wv();
            end
            if nargin > 1
                wv = bag.wv(:,k);
            else
                wv = bag.wv;
            end
        end

        % compute  image-word frequency
        function W = iwf(bag)

            N = bag.nimages;  % number of images

            % Create the word frequency matrix W
            %  column correspond to images
            %  row correspond to words
            %  each element is the number of occurences of that word in that iamge
            W = [];
            id = [bag.features.image_id];

            nl = bag.K - length(bag.stopwords);

            for i=1:bag.nimages
                % get the words associated with image i
                words = bag.words(id == i);

                % create columns of the W
                [w,f] = count_unique(words);
                v = zeros(nl,1);
                v(w) = f;
                W = [W v];
            end
        end

        function W = compute_wv(bag, bag2)

            if nargin == 2
                Wv = bag2.iwf();
                N = bag2.nimages;
                W = bag.iwf();
            else
                Wv = bag.iwf();
                N = bag.nimages;
                W = Wv;
            end

            Ni = sum( Wv'>0 );

            m = [];
            for i=1:bag.nimages
                % number of words in this image
                nd = sum( W(:,i) );

                % word occurrence frequency
                nid = W(:,i)';

                v = nid/nd .* log(N./Ni);
                v(~isfinite(v)) = 0;
                m = [m v'];
            end

            if nargout == 1
                W = m;
            else
                bag.wv = m;
            end
        end

        function [w,f] = wordfreq(bag)
        %BagOfWords.wordfreq Word frequency statistics
        %
        % [W,N] = B.wordfreq() is a vector of word labels W and the corresponding
        % elements of N are the number of occurrences of that word.
            [w,f] = count_unique(bag.words);
        end

        % compute similarity matrix
        function sim = similarity(bag1, bag2)
            wv1 = bag1.wordvector;
            wv2 = bag2.wordvector;
            for i=1:bag1.nimages
                for j=1:bag2.nimages
                    v1 = wv1(:,i); v2 = wv2(:,j);
                    sim(i,j) = dot(v1,v2) / (norm(v1) * norm(v2));
                end
            end
        end

        function display(bag)
        %BagOfWords.display Display value
        %
        % B.display() displays the parameters of the bag in a compact human
        % readable form.
        %
        % Notes::
        % - This method is invoked implicitly at the command line when the result
        %   of an expression is a BagOfWords object and the command has no trailing
        %   semicolon.
        %
        % See also BagOfWords.char.

            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            if loose
                disp(' ');
            end
            disp(char(bag))
            if loose
                disp(' ');
            end
        end

        function s = char(bag)
        %BagOfWords.char Convert to string
        %
        % S = B.char() is a compact string representation of a bag of words.
            s = sprintf(...
            'BagOfWords: %d features from %d images\n           %d words, %d stop words\n', ...
                length(bag.features), bag.nimages, ...
                bag.K-length(bag.stopwords), length(bag.stopwords));
        end

        function v = contains(bag, word)
        %BagOfWords.contains Find images containing word
        %
        % K = B.contains(W) is a vector of the indices of images in the sequence that
        % contain one or more instances of the word W.
            v = unique([bag.isword(word).image_id]);
        end

        function out = exemplars(bag, words, images, varargin)
        %BagOfWords.exemplars Display exemplars of words
        %
        % B.exemplars(W, IMAGES, OPTIONS) displays examples of the support regions of
        % the words specified by the vector W.  The examples are displayed as a table
        % of thumbnail images.  The original sequence of images from which the features
        % were extracted must be provided as IMAGES.
        %
        % IM = B.exemplars(W, IMAGES, OPTIONS) as above but returns the thumbnails
        % as a composite image.
        %
        % Options::
        % 'columns',N       Number of columns to display (default 10)
        % 'maxperimage',M   Maximum number of exemplars to display from any
        %                   one image (default 2)
        % 'width',W         Width of each thumbnail [pixels] (default 50)
        % 'label'           Display word labels on the thumbnails.

            opt.gap = 2;
            opt.columns = 10;
            opt.maxperimage = 2;
            opt.width = 50;
            opt.label = false;
            opt.rows = [];

            opt = tb_optparse(opt, varargin);

            % figure the number of exemplars to show, no more than opt.maxperimage
            % from any one image
            nexemplars = 0;
            exemplars = {};
            for w=words
                image_prev = [];
                count = 0;
                for f=bag.isword(w)
                    if f.image_id == image_prev
                        count = count + 1;
                        if count > opt.maxperimage
                            continue;
                        end
                    end

                    exemplars = [exemplars {{w, f}}];
                end
            end

            if isempty(opt.rows)
                nr = ceil( length(exemplars) / opt.columns);
            else
                nr = opt.rows;
            end
            nc = min(length(exemplars), opt.columns);

            n = min(length(exemplars), nr*nc);
            exemplars = exemplars(1:n);

            Ng = opt.width+opt.gap;
            composite = ones(nr*Ng, nc*Ng);

            % render the support regions into composite image
            row = 0; col = 0;
            for ex=exemplars
                ex = ex{1};
                word = ex{1}; f = ex{2};

                % extract it from the containing image
                support = f.support(images, opt.width);

                % paste it into the panel
                composite = ipaste(composite, support, [col row]*Ng, 'zero');

                % update row/column indices
                col = col + 1;
                if col >= opt.columns
                    row = row + 1;
                    col = 0;
                end
            end

            if nargout == 1
                % output specified, return the image
                out = composite;
            else
                if opt.label
                    % no output specified, optionally label the cells

                    if nargout == 0
                        idisp(composite, 'plain');
                    end

                    row = 0; col = 0;
                    for ex=exemplars
                        ex = ex{1};
                        word = ex{1}; f = ex{2};

                        text(col*Ng+opt.gap*2, row*Ng+3*opt.gap, ...
                            sprintf('%d #%d', word, f.image_id), 'Color', 'g')

                        % update row/column indices
                        col = col + 1;
                        if col >= opt.columns
                            row = row + 1;
                            col = 0;
                        end
                    end
                else
                    idisp(composite, 'plain');
                end
            end
        end
    end
end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function pf = loadPoleFigure_rigaku_txt2(fname,varargin)
% import data fom ana file
%
% Syntax
%   pf = loadPoleFigure_rigaku_txt2(fname)
%
% Input
%  fname  - filename
%
% Output
%  pf - @PoleFigure
%
% See also
% ImportPoleFigureData loadPoleFigure

try

  [data,~,~,~,header] = txt2mat(fname,'infoLevel',0);

  rho = data(:,1);
  step = rho(2)-rho(1);
  header = strsplit(header,'\n');
  id = find(~cellfun(@isempty,strfind(header,'Step')),1);
  step2 = cell2mat(textscan(header{id},'Step %d'));
  assert(step2(1) == step);

  data(:,1:2:end) = [];
  theta = ((size(data,2)-2):-1:0)*step;

  r = regularS2Grid('theta',theta*degree,'rho',rho*degree);

  assert(size(data,1)>5 && size(data,2)>5);

  pf = PoleFigure(string2Miller(fname),r,data(:,2:end),varargin{:});

catch
  interfaceError(fname);
end

end

author: mtex-toolbox

repo: mtex

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


text:
function binary = matlab_real_to_binary_string ( x )

%*****************************************************************************80
%
%% MATLAB_REAL_TO_BINARY_STRING converts a MATLAB real number to a binary string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the MATLAB real number.
%
%    Output, string BINARY, the binary string.
%

%
%  Convert X to a hexadecimal string.
%
  hex = matlab_real_to_hex_string ( x );
%
%  Convert hexadecimal to binary.
%
  binary = hex_to_binary_string ( hex );

  return
end


text:
function result  = gomp_spr(Y, K, L, epsilon, options)
    if nargin < 3
        L = 2;
    end
    if nargin < 4
        epsilon = 1e-3;
    end
    if nargin < 5
        options = struct;
    end
    if ~isfield(options, 'sparse_output')
        options.sparse_output = 1;
    end
    if ~isfield(options, 'verbose')
        options.verbose = 0;
    end
    result = mex_gomp_spr(Y, K, L, epsilon, ...
        options.sparse_output,...
        options.verbose);
end

author: indigits

repo: sparse-plex


text:
% Digital Video Stabilization and Rolling Shutter Correction using Gyroscopes
% Copyright (C) 2011 Alexandre Karpenko
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

function [dth, theta] = frame_rotations_non_lin(gyro, gyro_time, frame_time, t0, ts)
% computes the delta in theta between frame times, as well as the shear
% at frame time
dgt = diff(gyro_time);
theta = ((gyro(1:end-1,:) + gyro(2:end,:)) / 2) .* dgt(:,[1 1 1]);
theta = [0 0 0; cumsum(theta, 1)];

sigma2 = 4000;
gauss = exp(-(-120:120).^2 / sigma2);
gauss = gauss ./ sum(gauss);
gyro(:,1) = gyro(:,1) - conv(gyro(:,1), gauss, 'same');
gyro(:,2) = gyro(:,2) - conv(gyro(:,2), gauss, 'same');
gyro(:,3) = gyro(:,3) - conv(gyro(:,3), gauss, 'same');
th = ((gyro(1:end-1,:) + gyro(2:end,:)) / 2) .* dgt(:,[1 1 1]);
th = [0 0 0; cumsum(th, 1)];
dth = diff(lininterp(gyro_time + t0 + ts/2, th, frame_time));

author: alex-golts

repo: Video-Stabilization


text:
% erplab_interpolateElectrodes() - interpolate data channels
%
% Usage: EEG = erplab_interpolateElectrodes(inputEEG, replace_elecs, ignore_elecs, method)
%
% Inputs:
%     EEG           - EEGLAB dataset
%     replace_elecs - [integer array] indices of channels to interpolate.
%                     For instance, these channels might be bad.
%                     [chanlocs structure] channel location structure containing
%                     either locations of channels to interpolate or a full
%                     channel structure (missing channels in the current
%                     dataset are interpolated).
%     ignore_elecs  - Do not include these electrodes as input for interopolation
%     method        - [string] method used for interpolation (default is 'spherical').
%                     'invdist'/'v4' uses inverse distance on the scalp
%                     'spherical' uses superfast spherical interpolation.
% Output:
%     EEGOUT        - data set with bad electrode data replaced by
%                     interpolated data
%
% Author: Jason Arita
%
% Built off of EEGLAB's eeg_interp.m by Arnoud Delorme

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

function EEG = erplab_interpolateElectrodes(inputEEG, replace_elecs, ignore_elecs, method, varargin)

%% Error check input variables
if nargin < 2
    help erplab_interpolateElectrodes;
    return;
end;

% Error check: channel structure exists
tmplocs = inputEEG.chanlocs;
if isempty(tmplocs) || isempty([tmplocs.X])
    error(sprintf('Missing Channel Locations:\n\nChannel locations are needed to interpolate. Add channel locations through EEGLAB > Edit > Channel locations')); %#ok<*SPERR>
end

% Error check: Overlapping electrode lists in replace_elecss and
% ignore_elecs
overlap_elec = intersect(replace_elecs, ignore_elecs);
if ~isempty(overlap_elec)
    error(sprintf('Replace/Ignore Electrode Overlap:\n\nThere is overlap in the replace electrodes and the ignore electrodes.\nFix the input to avoid this overlap.\n\nBad electrodes:\t\t%s\nIgnore electrodes:\t%s\nOverlapping electrodes:\t%s', ...
        num2str(replace_elecs), ...
        num2str(ignore_elecs), ...
        num2str(overlap_elec)));
end

% Error check: When the user inputs electrodes not found in the input EEG dataset
missing_replace_elecss = setdiff(replace_elecs, [inputEEG.chanlocs.urchan]);
if(~isempty(missing_replace_elecss))
    error(sprintf('Missing Replace Electrodes:\n\nThe following electrodes entered as input were missing from the EEG dataset:\n\t\t%s', ...
        num2str(missing_replace_elecss)));
end

% Warning check: Missing ignored electrodes in EEG dataset
missing_ignore_elecs = setdiff(ignore_elecs, [inputEEG.chanlocs.urchan]);
if(~isempty(missing_ignore_elecs))
    warning(sprintf('Missing Ignored Electrodes:\n\nThe following electrodes entered as input were missing from the EEG dataset:\n\t%s', ...
        num2str(missing_ignore_elecs))); %#ok<SPWRN>
end

%%
EEG = inputEEG;

% If no method specified: Spherical interpolate default
if nargin < 4
    disp('Using spherical interpolation');
    method = 'spherical';
end

%% Handle optional input variables
optargs = {false}; % Default optional inputs

% Put defaults into the valuesToUse cell array,
% and overwrite the ones specified in varargin.
optargs(1:length(varargin)) = varargin;     % if vargin is empty, optargs keep their default values. If vargin is specified then it overwrites

% Place optional args into variable names
[displayEEG] = optargs{:};

%%
if isstruct(replace_elecs)

    % add missing channels in interpolation structure
    % -----------------------------------------------
    lab1 = { replace_elecs.labels };
    tmpchanlocs = EEG.chanlocs;
    lab2 = { tmpchanlocs.labels };
    [~, tmpchan] = setdiff_bc( lab2, lab1);
    tmpchan = sort(tmpchan);
    if ~isempty(tmpchan)
        newchanlocs = [];
        fields = fieldnames(replace_elecs);
        for index = 1:length(fields)
            if isfield(replace_elecs, fields{index})
                for cind = 1:length(tmpchan)
                    fieldval = getfield( EEG.chanlocs, { tmpchan(cind) },  fields{index});
                    newchanlocs = setfield(newchanlocs, { cind }, fields{index}, fieldval);
                end;
            end;
        end;
        newchanlocs(end+1:end+length(replace_elecs)) = replace_elecs;
        replace_elecs = newchanlocs;
    end;
    if length(EEG.chanlocs) == length(replace_elecs), return; end;

    lab1          = { replace_elecs.labels };
    tmpchanlocs   = EEG.chanlocs;
    lab2          = { tmpchanlocs.labels };
    [~, badchans] = setdiff_bc( lab1, lab2);

    fprintf('Interpolating %d channels...\n', length(badchans));

    if isempty(badchans), return; end;
    goodchans      = sort(setdiff(1:length(replace_elecs), badchans));

    % re-order good channels
    % ----------------------
    [~, tmp2, neworder] = intersect_bc( lab1, lab2 );
    [~, ordertmp2]      = sort(tmp2);
    neworder            = neworder(ordertmp2);
    EEG.data            = EEG.data(neworder, :, :);

    % looking at channels for ICA
    % ---------------------------
    %[~, sorti] = sort(neworder);
    %{ EEG.chanlocs(EEG.icachansind).labels; replace_elecs(goodchans(sorti(EEG.icachansind))).labels }

    % update EEG dataset (add blank channels)
    % ---------------------------------------
    if ~isempty(EEG.icasphere)

        [~, sorti] = sort(neworder);
        EEG.icachansind = sorti(EEG.icachansind);
        EEG.icachansind = goodchans(EEG.icachansind);
        EEG.chaninfo.icachansind = EEG.icachansind;

        % TESTING SORTING
        %icachansind = [ 3 4 5 7 8]
        %data = round(rand(8,10)*10)
        %neworder = shuffle(1:8)
        %data2 = data(neworder,:)
        %icachansind2 = sorti(icachansind)
        %data(icachansind,:)
        %data2(icachansind2,:)
    end;
    % { EEG.chanlocs(neworder).labels; replace_elecs(sort(goodchans)).labels }
    %tmpdata                  = zeros(length(replace_elecs), size(EEG.data,2), size(EEG.data,3));
    %tmpdata(goodchans, :, :) = EEG.data;

    % looking at the data
    % -------------------
    %tmp1 = mattocell(EEG.data(sorti,1));
    %tmp2 = mattocell(tmpdata(goodchans,1));
    %{ EEG.chanlocs.labels; replace_elecs(goodchans).labels; tmp1{:}; tmp2{:} }
    %EEG.data      = tmpdata;

    EEG.chanlocs  = replace_elecs;

else
    badchans  = replace_elecs;
    goodchans = setdiff_bc(1:EEG.nbchan, badchans);
    oldelocs  = EEG.chanlocs;
    EEG       = pop_select(EEG, 'nochannel', badchans);
    EEG.chanlocs = oldelocs;
    disp('Interpolating missing channels...');
end;

% find non-empty good channels
% ----------------------------
origoodchans  = goodchans;
chanlocs      = EEG.chanlocs;
nonemptychans = find(~cellfun('isempty', { chanlocs.theta }));

%% Remove ignored channels from interpolation computation
nonemptychans = setdiff_bc(nonemptychans, ignore_elecs);

[~, indgood ] = intersect_bc(goodchans, nonemptychans);
goodchans = goodchans( sort(indgood) );
datachans = getdatachans(goodchans,badchans);
badchans  = intersect_bc(badchans, nonemptychans);
if isempty(badchans), return; end;

%% Error check: Less than 2 good channels for interpolation
if length(goodchans) < 2
    feedback_str = [ ...
        sprintf('Cannot interpolate \t%.2d channels:\t', length(badchans)) ...
        sprintf(' %2.2d', badchans) ...
        '\n'...
        sprintf('using only\t\t%.2d channel(s):\t', length(goodchans)) ...
        sprintf(' %2.2d', goodchans) ...
        '\n' ...
        sprintf('skipping\t%.2d channels:\t', length(ignore_elecs))...
        sprintf(' %2.2d', ignore_elecs) ...
        '\n'];
    fprintf(feedback_str);
    error('Cannot interpolate with %d channels.\nEither make sure data set contains more than 2 channels \nAND that you are not ignoring too many channels', length(goodchans));
else
    % Print feedback to command window if error check passes
    feedback_str = [ ...
        sprintf('Interpolating \t%.2d channels:\t', length(badchans)) ...
        sprintf(' %2.2d', badchans) ...
        '\n'...
        sprintf('with\t\t%.2d channels:\t', length(goodchans)) ...
        sprintf(' %2.2d', goodchans) ...
        '\n' ...
        sprintf('skipping\t%.2d channels:\t', length(ignore_elecs))...
        sprintf(' %2.2d', ignore_elecs) ...
        '\n'];
    fprintf(feedback_str);
end

%% Warning check: Ignored bad channels
flatlinechans = intersect(replace_elecs, ignore_elecs);
if(~isempty(flatlinechans))
    warning('Warning: %d channels contain no data because they are interpolated with ignored channels: %s', length(flatlinechans), int2str(flatlinechans));
end

% scan data points
% ----------------
if strcmpi(method, 'spherical')
    % get theta, rad of electrodes
    % ----------------------------
    tmpgoodlocs = EEG.chanlocs(goodchans);
    xelec       = [ tmpgoodlocs.X ];
    yelec       = [ tmpgoodlocs.Y ];
    zelec       = [ tmpgoodlocs.Z ];
    rad         = sqrt(xelec.^2+yelec.^2+zelec.^2);
    xelec       = xelec./rad;
    yelec       = yelec./rad;
    zelec       = zelec./rad;
    tmpbadlocs  = EEG.chanlocs(badchans);
    xbad        = [ tmpbadlocs.X ];
    ybad        = [ tmpbadlocs.Y ];
    zbad        = [ tmpbadlocs.Z ];
    rad         = sqrt(xbad.^2+ybad.^2+zbad.^2);
    xbad        = xbad./rad;
    ybad        = ybad./rad;
    zbad        = zbad./rad;

    EEG.data    = reshape(EEG.data, EEG.nbchan, EEG.pnts*EEG.trials);
    %[~, tmp2, tmp3 tmpchans] = spheric_spline_old( xelec, yelec, zelec, EEG.data(goodchans,1));
    %max(tmpchans(:,1)), std(tmpchans(:,1)),
    %[~, tmp2, tmp3 EEG.data(badchans,:)] = spheric_spline( xelec, yelec, zelec, xbad, ybad, zbad, EEG.data(goodchans,:));
    [~, ~, ~, badchansdata] = spheric_spline( xelec, yelec, zelec, xbad, ybad, zbad, EEG.data(datachans,:));
    %max(EEG.data(goodchans,1)), std(EEG.data(goodchans,1))
    %max(EEG.data(badchans,1)), std(EEG.data(badchans,1))
    EEG.data = reshape(EEG.data, EEG.nbchan, EEG.pnts, EEG.trials);
elseif strcmpi(method, 'spacetime') % 3D interpolation, works but x10 times slower
    disp('Warning: if processing epoch data, epoch boundary are ignored...');
    disp('3-D interpolation, this can take a long (long) time...');

    tmpgoodlocs     = EEG.chanlocs(goodchans);
    tmpbadlocs      = EEG.chanlocs(badchans);
    [xbad ,ybad]    = pol2cart([tmpbadlocs.theta],[tmpbadlocs.radius]);
    [xgood,ygood]   = pol2cart([tmpgoodlocs.theta],[tmpgoodlocs.radius]);
    pnts            = size(EEG.data,2)*size(EEG.data,3);
    zgood           = 1:pnts;
    zgood           = repmat(zgood, [length(xgood) 1]);
    zgood           = reshape(zgood,prod(size(zgood)),1); %#ok<*PSIZE>
    xgood           = repmat(xgood, [1 pnts]);
    xgood           = reshape(xgood,prod(size(xgood)),1);
    ygood           = repmat(ygood, [1 pnts]);
    ygood           = reshape(ygood,prod(size(ygood)),1);
    tmpdata         = reshape(EEG.data, prod(size(EEG.data)),1);
    zbad            = 1:pnts;
    zbad            = repmat(zbad, [length(xbad) 1]);
    zbad            = reshape(zbad,prod(size(zbad)),1);
    xbad            = repmat(xbad, [1 pnts]); xbad = reshape(xbad,prod(size(xbad)),1);
    ybad            = repmat(ybad, [1 pnts]); ybad = reshape(ybad,prod(size(ybad)),1);
    badchansdata    = griddata(ygood, xgood, zgood, tmpdata,...
        ybad, xbad, zbad, 'nearest'); % interpolate data
else
    % get theta, rad of electrodes
    % ----------------------------
    tmpchanlocs     = EEG.chanlocs;
    [xbad ,ybad]    = pol2cart([tmpchanlocs( badchans).theta], ...
        [tmpchanlocs( badchans).radius]);
    [xgood,ygood]   = pol2cart([tmpchanlocs(goodchans).theta], ...
        [tmpchanlocs(goodchans).radius]);
    badchansdata = zeros(length(badchans), ...
        size(EEG.data,2)*size(EEG.data,3));
    fprintf('Points (/%d):', size(EEG.data,2)*size(EEG.data,3));

    for t=1:(size(EEG.data,2)*size(EEG.data,3)) % scan data points
        if mod(t,100) == 0, fprintf('%d ', t); end;
        if mod(t,1000) == 0, fprintf('\n'); end;

        %for c = 1:length(badchans)
        %   [h EEG.data(badchans(c),t)]= topoplot(EEG.data(goodchans,t),EEG.chanlocs(goodchans),'noplot', ...
        %        [EEG.chanlocs( badchans(c)).radius EEG.chanlocs( badchans(c)).theta]);
        %end;
        tmpdata = reshape(EEG.data, size(EEG.data,1), size(EEG.data,2)*size(EEG.data,3) );
        if strcmpi(method, 'invdist'), method = 'v4'; end;
        [~,~,badchansdata(:,t)] = griddata(ygood, xgood , double(tmpdata(datachans,t)'),...
            ybad, xbad, method); %#ok<GRIDD> % interpolate data
    end
    fprintf('\n');
end;

tmpdata               = zeros(length(replace_elecs), EEG.pnts, EEG.trials);
tmpdata(origoodchans, :,:) = EEG.data;
%if input data are epoched reshape badchansdata for Octave compatibility...
if length(size(tmpdata))==3
    badchansdata = reshape(badchansdata,length(badchans),size(tmpdata,2),size(tmpdata,3));
end
tmpdata(badchans,:,:) = badchansdata;
EEG.data = tmpdata;
EEG.nbchan = size(EEG.data,1);
EEG = eeg_checkset(EEG);

%% Display input EEG plot to user with rejection windows marked
% Needs to occur before actually deleting the time segments in EEG
if displayEEG
    % Plot EEG data with to-be-rejected time windows

    eegplotoptions = { ...
        'events',       EEG.event,          ...
        'srate',        EEG.srate,          ...
        'winlength',    10};

    % Display channel labels instead of numbers
    if ~isempty(EEG.chanlocs)
        eegplotoptions = [ eegplotoptions {'eloc_file', EEG.chanlocs}];
    end;

    eegplot(EEG.data, eegplotoptions{:});

end

%% Nested functions

function datachans = getdatachans(goodchans, badchans)
% Get data channels

datachans = goodchans;
badchans  = sort(badchans);
for index = length(badchans):-1:1
    datachans(find(datachans > badchans(index))) = ...
        datachans(find(datachans > badchans(index)))-1; %#ok<FNDSB>
end;

%% spherical splines
function [x, y, z, Res] = spheric_spline_old( xelec, yelec, zelec, values) %#ok<*DEFNU>

SPHERERES              = 20;
[x,y,z]                = sphere(SPHERERES);
x(1:(length(x)-1)/2,:) = []; x = x(:)';
y(1:(length(y)-1)/2,:) = []; y = y(:)';
z(1:(length(z)-1)/2,:) = []; z = z(:)';

Gelec = computeg(xelec,yelec,zelec,xelec,yelec,zelec);
Gsph  = computeg(x,y,z,xelec,yelec,zelec);

% equations are
% Gelec*C + C0  = Potential (C unknow)
% Sum(c_i) = 0
% so
%             [c_1]
%      *      [c_2]
%             [c_ ]
%    xelec    [c_n]
% [x x x x x]         [potential_1]
% [x x x x x]         [potential_ ]
% [x x x x x]       = [potential_ ]
% [x x x x x]         [potential_4]
% [1 1 1 1 1]         [0]

% compute solution for parameters C
% ---------------------------------
meanvalues = mean(values);
values = values - meanvalues; % make mean zero
C = pinv([Gelec;ones(1,length(Gelec))]) * [values(:);0];

% apply results
% -------------
Res = zeros(1,size(Gsph,1));
for j = 1:size(Gsph,1)
    Res(j) = sum(C .* Gsph(j,:)');
end
Res = Res + meanvalues;
Res = reshape(Res, length(x(:)),1);

function [xbad, ybad, zbad, allres] = spheric_spline( xelec, yelec, zelec, xbad, ybad, zbad, values)

newchans = length(xbad);
numpoints = size(values,2);

%SPHERERES = 20;
%[x,y,z] = sphere(SPHERERES);
%x(1:(length(x)-1)/2,:) = []; xbad = [ x(:)'];
%y(1:(length(x)-1)/2,:) = []; ybad = [ y(:)'];
%z(1:(length(x)-1)/2,:) = []; zbad = [ z(:)'];

Gelec = computeg(xelec,yelec,zelec,xelec,yelec,zelec);
Gsph  = computeg(xbad,ybad,zbad,xelec,yelec,zelec);

% compute solution for parameters C
% ---------------------------------
meanvalues = mean(values);
values = values - repmat(meanvalues, [size(values,1) 1]); % make mean zero

values = [values;zeros(1,numpoints)];
C = pinv([Gelec;ones(1,length(Gelec))]) * values;
clear values;
allres = zeros(newchans, numpoints);

% apply results
% -------------
for j = 1:size(Gsph,1)
    allres(j,:) = sum(C .* repmat(Gsph(j,:)', [1 size(C,2)]));
end
allres = allres + repmat(meanvalues, [size(allres,1) 1]);

% compute G function
% ------------------
function g = computeg(x,y,z,xelec,yelec,zelec)

unitmat = ones(length(x(:)),length(xelec));
EI = unitmat - sqrt((repmat(x(:),1,length(xelec)) - repmat(xelec,length(x(:)),1)).^2 +...
    (repmat(y(:),1,length(xelec)) - repmat(yelec,length(x(:)),1)).^2 +...
    (repmat(z(:),1,length(xelec)) - repmat(zelec,length(x(:)),1)).^2);

g = zeros(length(x(:)),length(xelec));
%dsafds
m = 4; % 3 is linear, 4 is best according to Perrin's curve
for n = 1:7
    if ismatlab
        L = legendre(n,EI);
    else % Octave legendre function cannot process 2-D matrices
        for icol = 1:size(EI,2)
            tmpL = legendre(n,EI(:,icol));
            if icol == 1, L = zeros([ size(tmpL) size(EI,2)]); end;
            L(:,:,icol) = tmpL;
        end;
    end;
    g = g + ((2*n+1)/(n^m*(n+1)^m))*squeeze(L(1,:,:));
end
g = g/(4*pi);

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%% Compare performance of gaimc to matlab_bgl
% While the |gaimc| library implements its graph routines in Matlab
% "m"-code, the |matlab_bgl| library uses graph algorithms from the Boost
% graph library in C++ through a mex interface.  Folklore has it that
% Matlab code with for loops like those required in the |gaimc| library is
% considerably slower.  This example examines this lore and shows that
% |gaimc| is typically within a factor of 2-4 of the mex code.

%%
% *This demo is unlikely to work on your own computer.*
% *They depend on having the MatlabBGL routines in one spot.*

%% Setup the environment
% We need MatlabBGL on the path
graphdir = '../graphs/';
matlabbgldir = '~/dev/matlab-bgl/4.0';
try
    addpath(matlabbgldir); % change this to your matlab_bgl path
    ci=components(sparse(ones(5)));
catch
    error('gaimc:performance_comparison','Matlab BGL is not working, halting...');
end

%%
% Check to make sure we are in the correct directory
cwd = pwd; dirtail = ['gaimc' filesep 'demo'];
if strcmp(cwd(end-length(dirtail)+1:end),dirtail) == 0
    error('%s should be executed from %s\n',mfilename,dirtail);
end

%%
% initalize the results structure
results=[];
mex_fast=0; mat_fast=0; mex_std=0; mat_std=0;

%% Depth first search
% To compare these functions, we have to make a copy and then delete it
copyfile(['..' filesep 'dfs.m'],'dfstest.m');
graphs = {'all_shortest_paths_example', 'clr-24-1', 'cs-stanford', ...
    'minnesota','tapir'};
nrep=30; ntests=100;
fprintf(repmat(' ',1,76));
for rep=1:nrep
    % Matlab needs 1 iteration to compile the function
    if nrep==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end
    for gi=1:length(graphs)
        load([graphdir graphs{gi} '.mat']); n=size(A,1);
        At=A'; [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai;
        for ti=1:ntests
            fprintf([repmat('\b',1,76) '%20s rep=%3i graph=%-30s trial=%4i'], ...
                'dfs',rep,graphs{gi},ti);
            v=ceil(n*rand(1));
            tic; d1=dfs(A,v); mex_std=mex_std+toc;
            tic; d2=dfs(At,v,struct('istrans',1,'nocheck',1));
              mex_fast=mex_fast+toc;
            tic; d3=dfstest(A,v); mat_std=mat_std+toc;
            tic; d4=dfstest(As,v); mat_fast=mat_fast+toc;
            if any(d1 ~= d2) || any(d2 ~= d3) || any(d3 ~= d4)
                error('gaimc:dfs','incorrect results from dijkstra');
            end
        end
    end
end
fprintf('\n');
delete('dfstest.m');
results(end+1).name='dfs';
results(end).mex_fast = mex_fast;
results(end).mat_fast = mat_fast;
results(end).mex_std = mex_std;
results(end).mat_std = mat_std;

%% Connected components
% To evaluate the performance of the connected components algorithm, we use
% sets of random graphs.
nrep=30;
szs=[1 10 100 5000 10000 50000];
comp_results=[mex_fast mat_fast mex_std mat_std];
for szi=1:length(szs)
    fprintf('\n%20s size=%-5i     ', 'scomponents', szs(szi));
    % Matlab needs 1 iteration to compile the function
    if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end
    for rep=1:nrep
        fprintf('\b\b\b\b'); fprintf(' %3i', rep);
        A=sprand(szs(szi),szs(szi),25/szs(szi));
        At=A'; [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai;
        tic; cc1=components(A); mex_std=mex_std+toc;
        tic; cc2=components(At,struct('istrans',1,'nocheck',1));
            mex_fast=mex_fast+toc;
        tic; cc3=scomponents(A); mat_std=mat_std+toc;
        tic; cc4=scomponents(As); mat_fast=mat_fast+toc;
        cs1=accumarray(cc1,1,[max(cc1) 1]);
        cs2=accumarray(cc2,1,[max(cc2) 1]);
        cs3=accumarray(cc3,1,[max(cc3) 1]);
        cs4=accumarray(cc4,1,[max(cc4) 1]);
        if any(cs1 ~= cs2) || any(cs2 ~= cs3) || any(cs2 ~= cs4)
            error('gaimc:scomponents','incorrect results from scomponents');
        end
    end
    comp_results(end+1,:) = [mex_fast mat_fast mex_std mat_std];
end
comp_results=diff(comp_results);
results(end+1).name='scomponents';
results(end).mex_fast = mex_fast;
results(end).mat_fast = mat_fast;
results(end).mex_std = mex_std;
results(end).mat_std = mat_std;

%%
% Plot the data for connected components
% This plot isn't too meaningful, so I've omitted it.

% plot(szs, comp_results,'.-');
% legend('mex fast','mat fast','mex std','mat std','Location','Northwest');
% ylabel('time'); xlabel('graph size');

%% Dijkstra's algorithm
% To evaluate the performance of Dijkstra's algorithm, we pick
graphs = {'clr-25-2', 'clr-24-1', 'cs-stanford', ...
    'minnesota', 'tapir'};
nrep=30; ntests=100; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0;
for rep=1:nrep
    for gi=1:length(graphs)
        load([graphdir graphs{gi} '.mat']); n=size(A,1);
        At=A'; [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai;
        for ti=1:ntests
            fprintf([repmat('\b',1,66) '%20s rep=%3i graph=%-20s trial=%4i'], ...
                'dijkstra',rep,graphs{gi},ti);
            v=ceil(n*rand(1));
            tic; d1=dijkstra_sp(A,v); mex_std=mex_std+toc;
            tic; d2=dijkstra_sp(At,v,struct('istrans',1,'nocheck',1));
              mex_fast=mex_fast+toc;
            tic; d3=dijkstra(A,v); mat_std=mat_std+toc;
            tic; d4=dijkstra(As,v); mat_fast=mat_fast+toc;
            if any(d1 ~= d2) || any(d2 ~= d3) || any(d3 ~= d4)
                error('gaimc:dijkstra','incorrect results from dijkstra');
            end
        end
    end
end
fprintf('\n');
results(end+1).name='dijkstra';
results(end).mex_fast = mex_fast;
results(end).mat_fast = mat_fast;
results(end).mex_std = mex_std;
results(end).mat_std = mat_std;

%% Directed Clustering coefficients
nrep=30; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0;
comp_results=[mex_fast mat_fast mex_std mat_std];
szs=[1 10 100 5000 10000 50000];
for szi=1:length(szs)
    fprintf('\n%20s size=%-5i     ', 'dirclustercoeffs', szs(szi));
    % Matlab needs 1 iteration to compile the function
    if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end
    for rep=1:nrep
        fprintf('\b\b\b\b'); fprintf(' %3i', rep);
        A=sprand(szs(szi),szs(szi),25/szs(szi));
        At=A';
        [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai;
        [cp ri ati]=sparse_to_csr(At); As.cp=cp; As.ri=ri; As.ati=ati;
        tic; cc1=clustering_coefficients(A); mex_std=mex_std+toc;
        tic; cc2=clustering_coefficients(At,struct('istrans',1,'nocheck',1));
            mex_fast=mex_fast+toc;
        tic; cc3=dirclustercoeffs(A); mat_std=mat_std+toc;
        tic; cc4=dirclustercoeffs(As); mat_fast=mat_fast+toc;
    end
    comp_results(end+1,:) = [mex_fast mat_fast mex_std mat_std];
end
fprintf('\n');
comp_results=diff(comp_results);
results(end+1).name='dirclustercoeffs';
results(end).mex_fast = mex_fast;
results(end).mat_fast = mat_fast;
results(end).mex_std = mex_std;
results(end).mat_std = mat_std;

%% Minimum spanning tree
nrep=30; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0;
comp_results=[];
szs=[10 100 5000 10000];
for szi=1:length(szs)
    fprintf('\n%20s size=%-5i     ', 'mst_prim', szs(szi));
    % Matlab needs 1 iteration to compile the function
    if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end
    for rep=1:nrep
        fprintf('\b\b\b\b'); fprintf(' %3i', rep);
        A=abs(sprandsym(szs(szi),25/szs(szi)));
        At=A';
        [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai;
        tic; T1=prim_mst(A); mex_std=mex_std+toc;
        tic; [t1i t1j t1v]=prim_mst(At,struct('istrans',1,'nocheck',1));
            mex_fast=mex_fast+toc;
        tic; T2=mst_prim(A,0); mat_std=mat_std+toc;
        tic; [t2i t2j t2v]=mst_prim(As,0); mat_fast=mat_fast+toc;
        T1f=sparse(t1i,t1j,t1v,size(A,1),size(A,2));
        T2f=sparse(t2i,t2j,t2v,size(A,1),size(A,2));
        if ~isequal(T1,T2) || ~isequal(T1f+T1f',T2f+T2f') || ...
                ~isequal(T1,T2f+T2f')
            keyboard
            warning('gaimc:mst_prim',...
                'incorrect results from mst_prim (%i,%i)', szi, rep);
        end
    end
    comp_results(end+1,:) = [mex_fast mat_fast mex_std mat_std];
end
fprintf('\n');
comp_results=diff(comp_results);
results(end+1).name='mst_prim';
results(end).mex_fast = mex_fast;
results(end).mat_fast = mat_fast;
results(end).mex_std = mex_std;
results(end).mat_std = mat_std;

%% Undirected Clustering coefficients
nrep=30; mex_fast=0; mat_fast=0; mex_std=0; mat_std=0;
undircc_results=[mex_fast mat_fast mex_std mat_std];
szs=[1 10 100 5000 10000 50000];
for szi=1:length(szs)
    fprintf('\n%20s size=%-5i     ', 'clustercoeffs', szs(szi));
    % Matlab needs 1 iteration to compile the function
    if szi==2, mex_fast=0; mat_fast=0; mex_std=0; mat_std=0; end
    for rep=1:nrep
        fprintf('\b\b\b\b'); fprintf(' %3i', rep);
        A=abs(sprandsym(szs(szi),25/szs(szi)));
        At=A';
        [rp ci ai]=sparse_to_csr(A); As.rp=rp; As.ci=ci; As.ai=ai;
        tic; cc1=clustering_coefficients(A,struct('undirected',1)); mex_std=mex_std+toc;
        tic; cc2=clustering_coefficients(At,struct('undirected',1,'istrans',1,'nocheck',1));
            mex_fast=mex_fast+toc;
        tic; cc3=clustercoeffs(A); mat_std=mat_std+toc;
        tic; cc4=clustercoeffs(As); mat_fast=mat_fast+toc;
    end
    undircc_results(end+1,:) = [mex_fast mat_fast mex_std mat_std];
end
%%
fprintf('\n');
undircc_results=diff(undircc_results);
results(end+1).name='clustercoeffs';
results(end).mex_fast = mex_fast;
results(end).mat_fast = mat_fast;
results(end).mex_std = mex_std;
results(end).mat_std = mat_std;

%% Summarize the results
% We are going to summarize the results in a bar plot based on the
% algorithm.  Each algorithm is a single bar
graphs = {'all_shortest_paths_example', 'clr-24-1', 'cs-stanford', ...
    'minnesota','tapir'};, where the performance of the
% mex code is 1.
nresults=length(results);
Ystd = zeros(nresults,1);
Yfast = zeros(nresults,1);
for i=1:nresults
    Ystd(i)=results(i).mat_std/results(i).mex_std;
    Yfast(i)=results(i).mat_fast/results(i).mex_fast;
end
bar(1:nresults,[Ystd Yfast]); set(gca,'XTickLabel',{results.name});
legend('Standard','Fast','Location','Northwest');

author: ckczzj

repo: CHAN

save_path: github-repos/MATLAB/ckczzj-CHAN


text:
function [vid_h_org, vid_l_bic] = load_video(dirFile,nFrames,upscale,part,param)

vid_h_org = [];
vid_l_bic = [];

for i = 1:nFrames
    filename = sprintf('original/Frame %03d.png', i);
    fullname = strcat(dirFile,filename);
    tmp = im2double(imread(fullname,'png'));
    [M N] = size(tmp(:,:,1));
    switch part
        case 1
            vid_h_org{i} = tmp;
        case 2
            vid_h_org{i} = imresize(tmp, 1/2);
        case 3 % half x half (288x352)
            vid_h_org{i} = tmp(1:M/2,1:N/2,:);
        case 4 % quarter x quarter (144x176)
            vid_h_org{i} = tmp(1:M/4,1:N/4,:);
        case 5 % 64x64
            vid_h_org{i} = tmp(65:128,1:64,:);
        case 6 % 32x64
            vid_h_org{i} = tmp(65:96,1:64,:);
    end
    vid_l_bic{i} = imresize(vid_h_org{i}, 1/upscale);
end

if param.SHOW_IMAGE
    figure(1), imshow(vid_h_org{1}); title('original (frame1)');
end

author: thomas-koehler

repo: SupER


text:
%% test1.m - exp03
% This script evaluates vector W (see test_bmrm.m and testValidatoin.m)
% against test data set (TST).
%
% 30-07-10 Michal Uricar
% 07-03-11 Michal Uricar, change to exp03

% clc;
close all; clearvars;

%% Timestamp

fprintf(1,'Started on %s\n\n', datestr(now));

%% Add path

addpath('./Functions/');
addpath('../../matlab_toolbox/mex/');

%% Select functions

load('./MAT/functions_gdisp.mat');

%% Load TST data

% load config
load('./MAT/config.mat');
lbpconfig_fname = ['./MAT/lbpconfig_nImages-' config.nImages '_hop-' num2str(config.heightOfPyramid) '.mat'];
images_fname = ['./MAT/ImagesTst_nImages-' config.nImages '.mat'];
gt_fname = ['./MAT/GTTst_nImages-' config.nImages '.mat'];

load(lbpconfig_fname);
load(images_fname);
load(gt_fname);
load('./results/options.mat');
S = options.S;

switch (config.nImages)
    case 'all'
        nImages = size(Images, 2);
    case 'half'
        nImages = size(Images, 2)/2;
    otherwise
        nImages = str2num(config.nImages);
end;

fprintf('Loading TST data...\n');
tic
data.lbp = lbpconfig; data.imSize = imSize; data.Images = Images; data.nImages = nImages;
options.S = S; data.Y = gt; data.options = options;
data.tmpData = getPsiMat(data, 1);
if (displacement)
    [data.options.PsiGS0, data.options.PsiGS1, data.options.PsiGS2] = getDisplacements(data.options);
end;
clear lbpconfig imSize Images nImages gt gt_fname images_fname lbpconfig_fname options S

% create mapTable for extracting components of W
data.mapTable = f_createMapTable(data.options, data.tmpData);
toc
fprintf('Done.\n');

w_fname = ['./results/W_nImages-' config.nImages '_hop-' num2str(config.heightOfPyramid) '.mat'];
load(w_fname);

load('./MAT/kappas.mat');

%% Test

% % start matlabpool
%try
%    matlabpool;
%catch ME1
%    idSegLast = regexp(ME1.identifier, '(?<=:)\w+$', 'match');
%    if (strcmp(idSegLast, 'OpenConnection'))
%        fprintf('Already connected to matlab pool. Skipping...\n');
%    end;
%end

fprintf('Running test...\n');

err_nose = inf(1, data.nImages);
err_canthus_rl = inf(1, data.nImages);
err_canthus_lr = inf(1, data.nImages);
err_mouth_corner_r = inf(1, data.nImages);
err_mouth_corner_l = inf(1, data.nImages);
err_canthus_rr = inf(1, data.nImages);
err_canthus_ll = inf(1, data.nImages);
err_mean = inf(1, data.nImages);
err_maxdist = inf(1, data.nImages);
L = zeros(1, data.nImages);

tic

for i = 1 : data.nImages
    [Test TestSparse] = getPsiMat(data, i);
    GT = data.Y{i};
    P = argmax_mex(data.options, W, data.mapTable, TestSparse);

    errors = sqrt(sum((P - GT).^2)) * kappa(i) * 100;

    err_nose(i) = errors(1);
    err_canthus_rl(i) = errors(2);
    err_canthus_lr(i) = errors(3);
    err_mouth_corner_r(i) = errors(4);
    err_mouth_corner_l(i) = errors(5);
    err_canthus_rr(i) = errors(6);
    err_canthus_ll(i) = errors(7);
    err_mean(i) = mean(errors);
    err_maxdist(i) = max(errors);

    L(i) = 1/data.options.M * sum( errors );
end;
Rtst = 1/data.nImages * sum(L);
toc

%try
%    matlabpool close;
%catch ME2
%    idSegLast = regexp(ME2.identifier, '(?<=:)\w+$', 'match');
%    if (strcmp(idSegLast, 'OpenConnection'))
%        fprintf('Cannot close matlabpool. Skipping...\n');
%    end;
%end;

%% Save errors

fprintf('Saving errors...\n');
options = data.options;
save('./results/errors.mat', 'options', 'L', 'err_nose', 'err_canthus_rl', ...
            'err_canthus_lr', 'err_mouth_corner_r', 'err_mouth_corner_l', ...
            'err_canthus_rr', 'err_canthus_ll', 'err_mean', 'err_maxdist', 'Rtst');
fprintf('Done.\n');

fprintf('Rtst = %.4f\n', Rtst);

%% Timestamp

fprintf(1,'Finished on %s\n\n', datestr(now));

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
% Butterworth IIR Filter
srate=1000; % sample rate
npts=2000; %npts in signal
Nyquist=srate/2; %Nyquist frequency
lpf=300; %low-pass frequency
order=5; %filter order
t=[0:npts-1]/srate; %time scale for plot
x=(rand(npts,1)*2)-1; % raw data from -1 to +1

[b,a]=butter(order,lpf/Nyquist); %create filter coefficients

filtered_data=filter(b,a,x); % filter using 'b' and 'a' coefficients

% Calculate FFT
fftpts=npts; % number points in FFT
hpts=fftpts/2; % half number of FFT points
binwidth=srate/fftpts;
f=[0:binwidth:srate-binwidth];
x_fft=abs(fft(x))/hpts; %scaled FFT of original signal
filtered_fft=abs(fft(filtered_data))/hpts; %scaled FFT of filtered signal

subplot(2,2,1)
plot(t,x);
title('Raw Time Series');
subplot(2,2,3)
plot(t,filtered_data);
title('Filtered Time Series');
xlabel('Time (s)');
subplot(2,2,2)
plot(f(1:hpts),x_fft(1:hpts));
title('Raw FFT');
subplot(2,2,4)
plot(f(1:hpts),filtered_fft(1:hpts));
title('Filtered FFT');
xlabel('Frequency (Hz)');

author: yandld


text:
function xscore = mergeAllClusters(rez)

Nfilt = size(rez.W,2);

xscore = zeros(Nfilt);
chckd  = zeros(Nfilt);

nmerges = 0;
for ik = 1:Nfilt
   s0 = rez.st3(:,2) ==ik;

   ipair = rez.iNeigh(:, ik);

   %%
   for j = 2:length(ipair)
       i2 = ipair(j);

       if chckd(ik, i2)
           continue;
       end

       j2 = find(rez.iNeigh(:, i2)==ik);

       if ~isempty(j2) && rez.simScore(ik, i2)>.5
          s1 = rez.st3(:,2) ==i2;

          c1 = rez.cProj(s0, [1 j]);
          c2 = cat(1, c1, rez.cProj(s1, [j2 1]));

          xs = c2(:,1) - c2(:,2);

          if ik==161 && i2==167
              keyboard;
          end

          xscore(i2, ik) = mergeScore(xs);
          chckd(i2, ik) = 1;

          if xscore(i2, ik)>0
              nmerges = nmerges + 1;
          end
       end
   end
   %%
   if rem(ik, 100)==1
       fprintf('Found %d merges, checked %d/%d clusters \n', nmerges, ik, Nfilt)
   end

end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function vanity
% VANITY    Display your File Exchange downloads and rank statistics,
%           add them to a history file (fex.mat in work directory, by
%           default), and play Handel's Hallelujah chorus if rank has
%           improved. Modify  code by entering URL of own author page
%           and location of history file, if desired.
% EXAMPLE : vanity (Oh, the FEX code metrics..)
% AUTHOR  : Dimitri Shvorob, dimitri.shvorob@vanderbilt.edu, 7/15/07
file = [matlabroot '\work\fex.mat'];
page = urlread('http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objectType=author&objectId=1095920');
d = findby('Downloads:\s(\d*)');
r = findby('Rank:\s(\d*)');
if ~exist(file,'file')
   h = [now d r];      %#ok
   save(file,'h')
   disp('History file created')
   fprintf('Current downloads/rank: %d/%d \n',d,r)
else
   load(file)
   fprintf('Former  downloads/rank: %d/%d \nCurrent downloads/rank: %d/%d \n',h(end,2),h(end,3),d,r)    %#ok
   if r < h(end,3)
      s = load('handel');
      sound(s.y,s.Fs)
   end
   h = [h; [now d r]]; %#ok
   save(file,'h')
end

function[n] = findby(pattern)
s = regexp(page,pattern,'tokens');
n = str2double(cell2mat(s{1}));
end

end

author: Sable

repo: mcbench-benchmarks


text:
function [XX,H] = projectedGrid ( P1, P2, P3, P4 , nx, ny);

% new formalism using homographies

a00 = [P1;1];
a10 = [P2;1];
a11 = [P3;1];
a01 = [P4;1];

% Compute the planart collineation:

[H] = compute_collineation (a00, a10, a11, a01);

% Build the grid using the planar collineation:

x_l = ((0:(nx-1))'*ones(1,ny))/(nx-1);
y_l = (ones(nx,1)*(0:(ny-1)))/(ny-1);

pts = [x_l(:) y_l(:) ones(nx*ny,1)]';

XX = H*pts;

XX = XX(1:2,:) ./ (ones(2,1)*XX(3,:));

repo: ekfmonoslam


text:
function C_hs = hdr_covariance(X,C,nh,nScans);
%
% C_hs = hdr_covariance(X,C,nh,nScans);
%
% Event-Related GLM tools:
% Compute the voxel-independent factor of the
% covariance matrix of a general linear model.
% This is taken from eq. (16) in the Greve
% theory paper (FS-FAST).
%
% original code by gb, 11/04:
% updated by ras, 05/05
C_hs = zeros(size(X,2));

for s = 1:nScans
    C_hs = C_hs + ((X(:,:,s)')*(C(:,:,s)^(-1))*X(:,:,s));
end

C_hs = C_hs^(-1);

return


text:
function test09 ( dim_num, n, z )

%*****************************************************************************80
%
%% TEST09 tests SPHERE_MEASURE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST09\n' );
  fprintf ( 1, '  SPHERE_MEASURE computes the SPHERE measure of quality.\n' );
  fprintf ( 1, '  Nonintersecting sphere volume    S = %14f\n', ...
    sphere_measure ( dim_num, n, z ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% eegplugin_4dneuroimaging() - EEGLAB plugin for importing 4D neuroimaging PDF files.
%
% Usage:
%   >> eegplugin_4dneuroimaging(fig, trystrs, catchstrs);
%
% Inputs:
%   fig        - [integer]  EEGLAB figure
%   trystrs    - [struct] "try" strings for menu callbacks.
%   catchstrs  - [struct] "catch" strings for menu callbacks.
%
% Create a plugin:
%   For more information on how to create an EEGLAB plugin see the
%   help message of eegplugin_besa() or visit http://www.sccn.ucsd.edu/eeglab/contrib.html
%
% Author: Christian Wienbruch, University of Konstanz, Clinical Psychology, 2005
%         (plugin interface: Arnaud Delorme, SCCN, INC, UCSD)

% Copyright (C) 2005 Christian Wienbruch, University of Konstanz, Clinical Psychology
%

function vers = eegplugin_4dneuroimaging(fig, trystrs, catchstrs)

  vers = '4dneuroimaging1.00';
  if nargin < 3
      error('eegplugin_4dneuroimaging requires 3 arguments');
  end;

  % add folder to path
  % ------------------
  if ~exist('pop_read4d')
      p = which('eegplugin_4dneuroimaging.m');
      p = p(1:findstr(p,'eegplugin_4dneuroimaging.m')-1);
      addpath([ p vers ] );
  end;

  % find import data menu
  % ---------------------
  menu = findobj(fig, 'tag', 'import data');

  % menu callbacks
  % --------------
  com = [ trystrs.no_check '[EEGTMP LASTCOM]= pop_read4d;' catchstrs.new_non_empty ];

  % create menus
  % ------------
  uimenu( menu, 'label',  'From 4D .m4d pdf file', 'callback', com, 'separator', 'on' );

author: PatternRecognition

repo: OpenBMI


text:
% eeg_point2lat() - convert latency in data points to latency in ms relative
%                   to the time locking. Used in eeglab().
% Usage:
%       >> [newlat ] = eeg_point2lat( lat_array, [], srate);
%       >> [newlat ] = eeg_point2lat( lat_array, epoch_array,...
%                                 srate, timelimits, timeunit);
% Inputs:
%   lat_array   - latency array in data points assuming concatenated
%                 data epochs (see eeglab() event structure)
%   epoch_array - epoch number corresponding to each latency value
%   srate       - data sampling rate in Hz
%   timelimits  - [min max] timelimits in 'timeunit' units (see below)
%   timeunit    - time unit in second. Default is 1 = seconds.
%
% Outputs:
%   newlat      - converted latency values (in 'timeunit' units) for each epoch
%
% Example:
%   tmpevent = EEG.event;
%   eeg_point2lat( [ tmpevent.latency ], [], EEG.srate, [EEG.xmin EEG.xmax]);
%   % returns the latency of all events in second for a continuous
%   % dataset EEG
%
%   eeg_point2lat( [ tmpevent.latency ], [ tmpevent.epoch ],
%                 EEG.srate, [EEG.xmin EEG.xmax]*1000, 1E-3);
%   % returns the latency of all events in millisecond for a dataset
%   % containing data epochs.
%
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2 Mai 2002
%
% See also: eeg_lat2point(), eeglab(), pop_editieventvals(), pop_loaddat()

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

function newlat = eeg_point2lat( lat_array, epoch_array, srate, timewin, timeunit);

if nargin <3
    help eeg_point2lat;
    return;
end;
if isempty( epoch_array )
    epoch_array = ones( size(lat_array) );
end;
if nargin <4
    timewin = 0;
end;
if nargin <5
	timeunit = 1;
end;

if length(lat_array) ~= length(epoch_array)
	if length(epoch_array)~= 1
		disp('eeg_point2lat: latency and epoch arrays must have the same length'); return;
	else
		epoch_array = ones(1,length(lat_array))*epoch_array;
	end;
end;
if length(timewin) ~= 2
    disp('eeg_point2lat: timelimits array must have length 2'); return;
end;
if iscell(epoch_array)
	epoch_array = [ epoch_array{:} ];
end;
if iscell(lat_array)
	lat_array = [ lat_array{:} ];
end

timewin = timewin*timeunit;

if length(timewin) == 2
    pnts = (timewin(2)-timewin(1))*srate+1;
else
    pnts = 0;
end;
newlat  = ((lat_array - (epoch_array-1)*pnts-1)/srate+timewin(1))/timeunit;
newlat = round(newlat*1E9)*1E-9;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip