text:
% function PairwiseMaskScreen(hfig,fMaskList_matchZ,fMask_matchZ_ZList);

% FC{i_fish,i_reg} = [];
% FC{i_fish,i_reg}.cIX = cIX;
% FC{i_fish,i_reg}.gIX = gIX;
% FC{i_fish,i_reg}.xyz = xyz;
% FC{i_fish,i_reg}.xyz_avr = xyz_avr;
% FC{i_fish,i_reg}.xyz_norm = xyz_norm;
% FC{i_fish,i_reg}.xyz_norm_avr = xyz_norm_avr;

%% Part 1: from zMasks.

range_fish = 1:18;
numFish = length(range_fish);

% screen ALL pairs within these candidate fMasks
% [instead of "for each zMask, find nodes"]

% for ref fish
%   load fMask_ref range;
%   for each ref mask
%     set fMask_ref;
%     for test fish,
%         load fMask_test range;
%         for each test mask,
%             load fMask_test;
%             compute perc overlap for this pair,pool;
% end;end;end;end

% treat ZBrain masks similarly to a ref fish

for i_fishnum = 1:numFish,% cycle through each fish as reference fish

    % choose reference fish
    i_fish_ref = range_fish(i_fishnum);

    % load fMask_ref range
    name = ['fMASKs_',num2str(i_fish),'.mat'];
    load(fullfile(data_masterdir,'fMASKs',name),'fMASKs');
    nfMasks = size(fMASKs.MaskDatabase,2);

end

%%
% function PairwiseClusterScreen
range_fish = 2:4;
numFish = length(range_fish);
%% Batch pre-load
numReg = 1;
i_reg = 1;
M_clus = cell(numFish,numReg);
for i_fish = range_fish,
    [cIX,gIX,M_xyz_norm] = GetDefaultClustersFromLoad(hfig,i_fish);
    M_clus{i_fish,i_reg} = [];
%     M_clus{i_fish,i_reg}.cIX = cIX;
    M_clus{i_fish,i_reg}.gIX = gIX;
    M_clus{i_fish,i_reg}.M_xyz_norm = M_xyz_norm;
end

%% screen ALL pairs (non-directional) between ALL fish
% init
TF_fishrange = cell(1,numFish);
Pairs_AllClusAllFish = []; % each row: {[fish#1, clusID1]}, {[fish#2, clusID2]}

% cycle through each fish as reference fish
for i_refnum = 1:numFish,
    i_fish_ref = range_fish(i_refnum);
    disp(['i_fish_ref = ' num2str(i_fish_ref)]);
    gIX_ref = M_clus{i_fish_ref,i_reg}.gIX;
    M_xyz_norm_ref = M_clus{i_fish_ref,i_reg}.M_xyz_norm;

    % init
    U_ref = unique(gIX_ref);
    numClus_ref = length(U_ref);
    TF_reffish = zeros(numClus_ref,numFish-1);

    % cycle through clusters in reference fish
    for i_clus_ref = 1:numClus_ref,
        if mod(i_clus_ref,10)==0,
            disp(['i_clus_ref = ' num2str(i_clus_ref)]);
        end
        clusID_ref = U_ref(i_clus_ref);
        IX = find(gIX_ref == clusID_ref);
        XYZ_ref = M_xyz_norm_ref(IX,:);

        % cycle through test fish for given ref fish
        for i_testnum = (i_refnum+1):numFish,% only do ordered pairwise test (other half redundant)
            i_fish_test = range_fish(i_testnum);
            gIX_test = M_clus{i_fish_test,i_reg}.gIX;
            M_xyz_norm_test = M_clus{i_fish_test,i_reg}.M_xyz_norm;

            U_test = unique(gIX_test);
            numClus_test = length(U_test);
            scores_test = zeros(numClus_test,1);

            % cycle through clusters in test fish
            for i_clus_test = 1:numClus_test,
                clusID_test = U_test(i_clus_test);
                IX = find(gIX_test == clusID_test);
                XYZ_test = M_xyz_norm_test(IX,:);

                % compute distance score
                scores_test(i_clus_test,1) = ClusterDistanceTestD12(XYZ_ref,XYZ_test);
            end

            % save qualifying pairs for this test fish
            IX_pass = find(scores_test);
            TF_reffish(i_clus_ref,i_testnum) = length(IX_pass);
            if ~isempty(IX_pass),
                clusID_test_pass = U_test(IX_pass);
                for i = 1:length(clusID_test_pass),
                    Pairs_AllClusAllFish = [Pairs_AllClusAllFish; ...
                        {[i_fish_ref, clusID_ref]}, {[i_fish_test, clusID_test_pass(i)]}]; %#ok<AGROW>
                end
            end
        end
    end
    TF_fishrange{i_refnum} = TF_reffish;
end

author: xiuyechen

repo: FishExplorer

path: github-repos/MATLAB/xiuyechen-FishExplorer/FishExplorer-c61392cf0835480d64fc03c15f1992935fdc7106/unused analyses/conserved cluster screens/obsl mask screen/PairwiseMaskScreen.m


text:
%% Copyright (C) 2014-2016, 2019-2020, 2022 Colin B. Macdonald
%% Copyright (C) 2017 Mike Miller
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @deftypefun {[@var{A}, @var{db}] =} python_ipc_driver (@dots{})
%% Private function: run Python/SymPy command and return objects.
%%
%% @var{A} is the resulting object, which might be an error code.
%%
%% @var{db} usually contains diagnostics to help with debugging
%% or error reporting.
%% @end deftypefun

function [A, db] = python_ipc_driver(what, cmd, varargin)

  which_ipc = sympref('ipc');

  %% version check
  if exist('OCTAVE_VERSION', 'builtin')
    if (compare_versions (OCTAVE_VERSION (), '5.1', '<'))
      fprintf(['\n********************************************************************\n' ...
               'Your Octave version is %s but Octsympy is currently untested on\n' ...
               'anything before 5.1.\n\n'], OCTAVE_VERSION ())
      warning('Old Octave version detected: probably trouble ahead!')
      fprintf('\n********************************************************************\n\n')
      fflush(stdout);
    end
  end

  if (strcmp(lower(which_ipc), 'default'))
    % TODO: may need to adjust, ideally just 'py'
    if (exist ('pyversion') && exist ('pyexec') && exist ('pyeval'))
      which_ipc = 'native';
    elseif (exist ('popen2') > 1)
      which_ipc = 'popen2';
    else
      which_ipc = 'system';
    end
  end

  switch lower(which_ipc)
    case 'native'
      [A, db] = python_ipc_native(what, cmd, varargin{:});

    case 'popen2'
      [A, db] = python_ipc_popen2(what, cmd, varargin{:});

    case 'system'
      if (ispc () && (~isunix ()))
        [A, db] = python_ipc_sysoneline(what, cmd, varargin{:});
      else
        [A, db] = python_ipc_system(what, cmd, false, varargin{:});
      end

    case 'systmpfile'
      %% for debugging, not intended for long-term usage
      [A, db] = python_ipc_system(what, cmd, true, varargin{:});

    case 'sysoneline'
      %% for debugging, not intended for long-term usage
      [A, db] = python_ipc_sysoneline(what, cmd, varargin{:});

    otherwise
      if (strcmp (what, 'reset'))
        A = true;
        db = [];
      else
        error ('invalid ipc mechanism')
      end
  end

repo: octsympy


text:
function h = histc2( A, edges, wtMask )
% Multidimensional histogram count with weighted values.
%
% Creates a histogram h of the values in A [n x nd], with edges as
% specified. If nd==1, that is when A is a vector, the resulting histogram
% h is a vector of length nBins, where nBins=length(edges)-1. h(q) contains
% the weighted count of values v in A such that edges(q) <= v < edges(q+1).
% h(nBins) additionally contains the weighted count of values in A such
% that v==edges(nBins+1) -- which is different then how histc treates the
% boundary condition. Finally, h is normalized so that sum(h(:))==1.
%
% It usually makes sense to specify edges explicitly, especially if
% different histograms are going to be compared.  In general, edges must
% have monotonically non-decreasing values.  Also, if the exact bounds are
% unknown then it is convenient to set the first element in edges to -inf
% and the last to inf.  If h = histc2( A, nBins, ...), edges are
% automatically generated and have bins equally spaced between min(A) and
% max(A). That is the edges vector is generated by:
%  edges = linspace( min(A)-eps, max(A)+eps, nBins+1 );
%
% If nd>1, that is when A is a 2d matrix instead of a vector, the created
% histogram is multi-dimensional with dimensions nBins^nd, where each bin
% h(q1,...,qnd) contains the the weighted count of vectors v in A such that
% edges{k}(qk) <= v(k) < edges{k}(qk+1), for k=1,...,nd.  Note that if nd>1
% edges may be a cell vector where each element is a vector of edges or a
% scalar nBins as before.
%
% Each value in A may optionally have an associated weight given by wtMask,
% which should have the same number of elements as A. If not specified, the
% default is wtMask=ones(n,1).
%
% USAGE
%  h = histc2( A, edges, [wtMask] )
%
% INPUTS
%  A           - [n x nd] 2D numeric array
%  edges       - quantization bounds, see above
%  wtMask      - [] length [n] vector of weights
%
% OUTPUTS
%  h           - nd histogram [nBins^nd]
%
% EXAMPLE - 1D histograms
%  A=filterGauss([1000 1000],[],[],0); A=A(:); n=length(A);
%  h1 = histc2( A, 25 );              figure(1); bar(h1);
%  h2 = histc2( A, 25, ones(n,1) );   figure(2); bar(h2);
%  h3 = histc2( A, 25, A );           figure(3); bar(h3);
%
% EXAMPLE - 2D histograms
%  A=filterGauss([1000 1000],[],[],0); A=A(:); n=length(A);
%  h=histc2( [A A], 25 );    figure(1); im(h);  % decreasing along diag
%  h=histc2( [A A], 25, A ); figure(2); im(h);  % constant along diag
%
% See also HISTC, ASSIGNTOBINS, BAR
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 ); wtMask=[]; end;
if( ~isa(A,'double') ); A=double(A); end;
if( ~ismatrix(A) ); error('A must be a 2 dim array'); end;
[n,nd] = size(A);
if( ~isempty(wtMask) && n~=numel(wtMask) )
  error( 'wtMask must have n elements (A is nxnd)' ); end

if( nd==1 )
  % if nBins given instead of edges calculate edges
  if(length(edges)==1)
    edges = linspace(min(A)-eps,max(A)+eps,edges+1);
  end

  % create 1d histogram
  if(isempty(wtMask))
    h = histc( A, edges );
    h(end-1) = h(end-1)+h(end);
    h = h(1:end-1); h = h / sum(h);
  else
    h = histc2c( A, wtMask, edges );
    h = h / sum(h);
  end

else
  % if nBins given instead of edges calculate edges per dimension
  if( ~iscell(edges ) )
    edges=repmat({edges},[1 nd]);
  elseif( length(edges)~=nd )
    error( 'Illegal dimensions for edges' );
  end
  for i=1:length( edges );
    if(length(edges{i})==1)
      edges{i}=linspace(min(A(:,i))-eps,max(A(:,i))+eps,edges{i}+1);
    end
  end

  % create multidimensional histogram
  if( isempty(wtMask) ); wtMask=ones(1,n); end;
  h = histc2c( A, wtMask, edges{:} );
  h = h / sum(h(:));
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
%% Snake Demo
%
% An example using a SNAKETOY function.
%
% This function creates a snake toy graphic.  Each block
% in the snake is an hgtransform object.  A vector of angles (in
% radians) can be passed in to create a new shape by twisting the
% snake.  Twisting occurs using the hgtransform's 'Matrix'
% property.
%
% <snaketoy.m>

% Copyright (C) 2005-2011 The MathWorks Inc.

% Create a snake
snake = snaketoy;

%% Transform the snake
%
% By passing in a list of angles, you can transform the snake into
% a new shape.
snaketoy(snake,[ 0 0 0 pi pi 0 0 0 pi 0 pi pi ...
                   0 pi 0 pi pi 0 pi 0 pi pi 0 pi]);
camzoom(2)
%% Named Shapes
%
% Some shapes have pre-defined names encoded in the file.
% You can reference these by name.
snaketoy(snake, 'box');

%% Transformation Methods
%
% You can cause a shape to transform in different ways by moving
% linearly through the pattern from one of the snake to the other,
% choosing random joints to twist, or twisting all joints at the
% same time.
snaketoy(snake, 'dog','once');

%% Demonstration Mode
%
% You can run the script in demo mode by calling SNAKETOY
% with no arguments.  This feature runs infinitly, constantly
% transforming the snake into new shapes using arbitrary
% transformation methods.

% snaketoy

author: Sable

repo: mcbench-benchmarks


text:
function y = cat(varargin)
%CAT (overloaded)

switch varargin{1}
    case 1
        y = vertcat(varargin{2:end});
    case 2
        y = horzcat(varargin{2:end});
    otherwise
        for i = 2:nargin
            if isa(varargin{i},'sdpvar');
                varargin{i} = ndsdpvar(varargin{i});
            end
        end
        y = cat(varargin{:});
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function []=make_html_contents
disp('note that we assume to be in a directoy called spider');
'going..'

 ss=spider_path;

f=fopen('core_algorithms.txt','rt');
fseek(f,0,'eof');
l=ftell(f);
fseek(f,0,'bof');
core_algos=char(fread(f,l)');
fclose(f);

 fin = fopen([ss 'Contents.m'],'r')
 F = fread(fin);
 s = char(F');
 fclose(fin);

% fid = fopen(['c:\spider/doc/web/objects.html'],'w')
 fid = fopen(['objects.html'],'w')
 fprintf(fid,'\n');
 fprintf(fid,'\n');
 fprintf(fid,'\n');
 a='\n'; fprintf(fid,a);
a='\n'; fprintf(fid,a);
a='\n'; fprintf(fid,a);
a='\n'; fprintf(fid,a);
a='\n'; fprintf(fid,a);
a='
 \n'; fprintf(fid,a);
a=' \n'; fprintf(fid,a);
a='<title> SPIDER </title>\n'; fprintf(fid,a);
a=' \n'; fprintf(fid,a);
a='\n'; fprintf(fid,a);
a=' The Spider Objects\n'; fprintf(fid,a);
a=' \n'; fprintf(fid,a);
a='\n'; fprintf(fid,a);
a='
\n'; fprintf(fid,a);
a=' \n'; fprintf(fid,a);
a='
\n'; fprintf(fid,a);
a='
\n'; fprintf(fid,a);
a='
\n'; fprintf(fid,a);
a='\n'; fprintf(fid,a);
a=' \n'; fprintf(fid,a);
a='\n'; fprintf(fid,a);
a='  \n'; fprintf(fid,a);
a='
\n'; fprintf(fid,a);

f=[find(s=='%')]; s(f)=' '; % purge comments
s=s(find(s~=10));
f=[find(double(s)==13)]; % find line feeds
for i=2:length(f)-1
   process(s(f(i):f(i+1)-1),fid,core_algos);
end
a='
\n'; fprintf(fid,a);
fprintf(fid,' Objects marked with - are available in the extra package. ')
fclose(fid);

function a=process(s,fid,core_algos)
%------------------------
a=[s(1:length(s)) ]; doit=1;

%a; keyboard

%%------ search for title ----------
if findstr('objects',a) & not(s(4)==' ')
    a=['

  	    ' char(s) '

'];
    fprintf(fid,'\n');
    fprintf(fid,'%s',a);
    doit=0;
end
%%----------------------------------

if doit  %% -----add font colors
    f=min(find(a=='-'));
    if not(isempty(f))
        a=s(2:f-1); a2=s(f+1:length(s));
        a=strtrim(a);
        a2=strtrim(a2);
        f1=min(find(not(a==' ')));
        f2=max(find(not(a==' ')));
        %f1=min(find((a>='a' & a<='z')));
        %f2=max(find((a>='a' & a<='z')));
        a_bef=a(1:f1-1); a_aft=a(f2+1:length(a));
        a=a(f1:f2);
        obj_a=a;

%         make_html_help(a);
        a=[a_bef '' a ' ' a_aft ];
        s=['

  ' a '	 ' ispartofcore(obj_a,core_algos) '	  ' a2 char(13) '

'];
        fprintf(fid,'%s',s);
    end
%    fprintf(fid,'
\n');
end

a=1; return;

function s=ispartofcore(a,core_algos)

if(isempty(strfind(core_algos,a)))
    s='-';
else
    s='+';
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [CL] = visutil_getCommonRange(erp,ivals,varargin)
%VISUTIL_GETMAXRANGE - get maximum CLim range for several scalp plots of
%data in erp defined by ivals
%
%Synopsis:
% CLIM= visutil_getMaxRange(erp, ivals)
%
%Input:
% erp: eeg data
% ERP: struct of epoched EEG data.
% IVAL: time intervals for which scalp topography are to be plotted.

% irene, 10/2015
props= {
        'Class',     [],   '';
        'CLim',                 'sym',             'CHAR|DOUBLE[2]'};

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

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
nIvals=size(ivals,1);

for ii= 1:nIvals,
    if any(isnan(ivals(ii,:))),
      continue;
    end
    eee= erp;
    if nargin>=2 && ~isempty(ivals(ii,:)) && ~sum(any(isnan(ivals))),
         eee= proc_selectIval(eee, ivals(ii,:), 'IvalPolicy','minimal');
    end
    if ~isempty(opt.Class),
         eee= proc_selectClasses(eee, opt.Class);
    end
    if max(sum(eee.y,2))>1,
        eee= proc_average(eee);
    end
    eee.x= mean(eee.x,1);
    min_max(ii,1)=min(min(eee.x));
    min_max(ii,2)=max(max(eee.x));
end

if isequal(opt.CLim, 'sym'),
  zgMax= max(abs(min_max(:,2)));
  CL= [-zgMax zgMax];
elseif isequal(opt.CLim, 'range'),
  CL= [min(min_max(:,1)) max(min_max(:,2))];
elseif isequal(opt.CLim, '0tomax'),
  CL= [0.0001*diff([min(min_max(:,1)) max(min_max(:,2))]) max(min_max(:,2))];
elseif isequal(opt.CLim, 'minto0'),
  CL= [min(min_max(:,1)) 0.0001*diff([min(min_max(:,1)) max(min_max(:,2))])];
elseif isequal(opt.CLim, 'zerotomax'),
  CL= [0 max(min_max(:,2))];
elseif isequal(opt.CLim, 'mintozero'),
  CL= [min(min_max(:,1)) 0];
end
if diff(CL)==0, CL(2)= CL(2)+eps; end
end

author: bbci


text:
function [ pairwise ] = load_pairwise_data( p )

pairwise = struct();
pairwise.graph = [];
nextreg = p.nextreg;

if nextreg
    load(p.pairwise_relations, 'means', 'std_devs', 'graph');
    std_devs = permute(std_devs, [3 4 1 2]);
    means = permute(means, [3 4 1 2]);
    pairwise.graph = graph;
    pairwise.means = means;
    pairwise.std_devs = std_devs;
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function b = sspsl ( ap, n, kpvt, b )

%*****************************************************************************80
%
%% SSPSL solves the real symmetric system factored by SSPFA.
%
%  Discussion:
%
%    To compute inverse(A) * C where C is a matrix with P columns:
%
%      call sspfa ( ap, n, kpvt, info )
%
%      if ( info /= 0 ) go to ...
%
%      do j = 1, p
%        call sspsl ( ap, n, kpvt, c(1,j) )
%      end do
%
%    A division by zero may occur if SSPCO has set RCOND == 0.0D+00
%    or SSPFA has set INFO /= 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real ( kind = 8 ) AP(N*(N+1)/2), the output from SSPFA.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer KPVT(N), the pivot vector from SSPFA.
%
%    Input, real B(N), the right hand side.
%
%    Output, real B(N), the solution.
%

%
%  Loop backward applying the transformations and D inverse to B.
%
  k = n;
  ik = floor ( ( n * ( n - 1 ) ) / 2 );

  while ( 0 < k )

    kk = ik + k;

    if ( 0 <= kpvt(k) )
%
%  1 x 1 pivot block.
%
      if ( k ~= 1 )

        kp = kpvt(k);
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end
%
%  Apply the transformation.
%
        b(1:k-1) = saxpy ( k-1, b(k), ap(ik+1:ik+k-1), 1, b(1:k-1), 1 );

      end
%
%  Apply D inverse.
%
      b(k) = b(k) / ap(kk);
      k = k - 1;
      ik = ik - k;

    else
%
%  2 x 2 pivot block.
%
      ikm1 = ik - ( k - 1 );

      if ( k ~= 2 )

        kp = abs ( kpvt(k) );
%
%  Interchange.
%
        if ( kp ~= k-1 )
          temp = b(k-1);
          b(k-1) = b(kp);
          b(kp) = temp;
        end
%
%  Apply the transformation.
%
        b(1:k-2) = saxpy ( k-2, b(k), ap(ik+1:ik+k-2), 1, b(1:k-2), 1 );
        b(1:k-2) = saxpy ( k-2, b(k-1), ap(ikm1+1:ikm1+k-2), 1, b(1:k-2), 1 );

      end
%
%  Apply D inverse.
%
      km1k = ik + k - 1;
      kk = ik + k;
      ak = ap(kk) / ap(km1k);
      km1km1 = ikm1 + k - 1;
      akm1 = ap(km1km1) / ap(km1k);
      bk = b(k) / ap(km1k);
      bkm1 = b(k-1) / ap(km1k);
      denom = ak * akm1 - 1.0;
      b(k) = ( akm1 * bk - bkm1 ) / denom;
      b(k-1) = ( ak * bkm1 - bk ) / denom;
      k = k - 2;
      ik = ik - ( k + 1 ) - k;

    end

  end
%
%  Loop forward applying the transformations.
%
  k = 1;
  ik = 0;

  while ( k <= n )

    if ( 0 <= kpvt(k) )
%
%  1 x 1 pivot block.
%
      if ( k ~= 1 )
%
%  Apply the transformation.
%
        b(k) = b(k) + sdot ( k-1, ap(ik+1:ik+k-1), 1, b(1:k-1), 1 );
        kp = kpvt(k);
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end

      end

      ik = ik + k;
      k = k + 1;

    else
%
%  2 x 2 pivot block.
%
      if ( k ~= 1 )
%
%  Apply the transformation.
%
        b(k) = b(k) + sdot ( k-1, ap(ik+1:ik+k-1), 1, b(1:k-1), 1 );
        ikp1 = ik + k;
        b(k+1) = b(k+1) + sdot ( k-1, ap(ikp1+1:ikp1+k-1), 1, b(1:k-1), 1 );
        kp = abs ( kpvt(k) );
%
%  Interchange.
%
        if ( kp ~= k )
          temp = b(k);
          b(k) = b(kp);
          b(kp) = temp;
        end

      end

      ik = ik + k + k + 1;
      k = k + 2;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = acot(a)
%ACOT         Gradient inverse cotangent acot(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 04/07/04     S.M. Rump  accelaration for sparse input
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  wng = warning;
  warning off

  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  ax = (-1) ./ ( 1 + sqr(full(a.x(:))) );
  a.x = acot(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  warning(wng);

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function H = recnotch(notch,mode,M,N,W,S)
%RECNOTCH Generates axes notch filter transfer functions.
%   H = RECNOTCH(NOTCH,MODE,M,N,W,S) generates an M-by-N notch filter
%   transfer function consisting of symmetric pairs of rectangles of
%   width W placed on the vertical and/or horizontal axes of the
%   (centered) frequency rectangle. W must be an odd integer
%   to preserve the symmetry of the filtered Fourier transform.
%
%   NOTCH can be:
%
%      'reject'    Notchreject filter transfer function.
%
%      'pass'      Notchpass filter transfer function.
%
%   MODE can be:
%
%      'vertical'     Filtering on vertical axis only.
%
%      'horizontal'   Filtering on horizontal axis only.
%
%      'both'         Filtering on both axes.
%
%   If MODE is 'vertical', then vertical rectangles start at +/- S
%   pixels from the center of the frequency rectangle along the vertical
%   axis and extend to both ends of that axis. Similarly, if MODE is
%   'horizontal', then horizontal rectangles start at +/- S pixels from
%   the center and extend to both ends of that axis. If MODE is 'both',
%   then rectangles are placed in both directions starting at +/- S.
%   Alternatively, if MODE is 'both', then S can be a two-element
%   vector, [SV SH], specifying the starting rectangle locations in both
%   directions.
%
%   H = RECNOTCH(NOTCH,MODE,M,N) uses W = 3 and S = 1.
%
%   H is of floating point class double. It is returned uncentered for
%   consistency with the filtering function dftfilt. To view H as an
%   image or mesh plot, it should be centered using Hc = fftshift(H).
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Defaults.
if nargin < 6
   S = 1;
end
if nargin < 5
   W = 3;
end
if ~isodd(W)
   error('W must be odd.')
end

% AV and AH are rectangle amplitude values for the vertical and
% horizontal rectangles: 0 for notchreject and 1 for notchpass. The
% transfer functions are computed initially as reject transfer functions
% and then changed to pass if so specified in NOTCH.
if isequal(mode,'vertical')
   % Reject only along vertical axis.
   AV = 0;
   AH = 1;
   SV = S;
   SH = 1;
elseif isequal(mode,'horizontal')
   % Reject only along horizontal axis.
   AV = 1;
   AH = 0;
   SH = S;
   SV = 1;
elseif isequal(mode,'both')
   % Reject along both axes.
   AV = 0;
   AH = 0;
   if numel(S) == 1
      SV = S;
      SH = S;
   else
      SV = S(1);
      SH = S(2);
   end
else
   error('Unknown mode')
end

% Begin computation. The transfer function is generated as a 'reject'
% function. At the end, it is changed to a 'pass' function if it is so
% specified in parameter NOTCH.
H = rectangleReject(M,N,W,SV,SH,AV,AH);

% Finished computing the rectangle notch transfer function. Format the
% output so it is returned uncentered for use in function dftfilt, and
% change H to a 'pass' function if so specified.
H = ifftshift(H);
if isequal(notch,'pass')
   H = 1 - H;
end

%----------------------------------------------------------------------%
function H = rectangleReject(M,N,W,SV,SH,AV,AH)
% Preliminaries.
H = ones(M,N);
% Center of frequency rectangle.
UC = floor(M/2) + 1;
VC = floor(N/2) + 1;
% Width limits.
WL = (W - 1)/2;
% Compute rectangle notches with respect to center.
% Left, horizontal rectangle.
H(UC-WL:UC+WL, 1:VC-SH) = AH;
% Right, horizontal rectangle.
H(UC-WL:UC+WL, VC+SH:N) = AH;
% Top vertical rectangle.
H(1:UC-SV, VC-WL:VC+WL) = AV;
% Bottom vertical rectangle.
H(UC+SV:M, VC-WL:VC+WL) = AV;

author: dipum

repo: dipum-toolbox

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


text:
function [] = plotIntrinsicStatsAll(cameraCBparametersAllCams)
%% function for plotting intrinsics statistics for all cameras in STEP0
%
% INPUTS:
% * cameraCBparametersAllCams: a cell array containing cameraCBparameters structures for all cameras
%
% OUTPUTS:
% * a figure plotting intrinsics statistics for all cameras
% tab 1-  plot original parameters
% tab 2- plot parameters after distortion
% tab 3- original and corrected together

%%
hf=figure('units','normalized','outerposition',[.1 .1 .8 .8]);

Ncam=numel(cameraCBparametersAllCams);

tabgp = uitabgroup(hf);

%% BEFORE CORRECTION

tab(1) = uitab(tabgp,'Title','Before distortion correction');
axes('Parent',tab(1));

for ic=1:Ncam
    RadialDistortion(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.RadialDistortion;
    RadialDistortionError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.RadialDistortionError;
    TangentialDistortion(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.TangentialDistortion;
    TangentialDistortionError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.TangentialDistortionError;
    Skew(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.Skew;
    SkewError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.SkewError;
    FocalLength(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.FocalLength;
    FocalLengthError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.FocalLengthError;
    PrincipalPoint(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.PrincipalPoint;
    PrincipalPointError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.PrincipalPointError;
    imageSizeAll(ic,:)=cameraCBparametersAllCams{ic}.imagesInfo.imageSize;
    camInds(ic)=cameraCBparametersAllCams{ic}.icam;
end

hf.Name=['camera intrinsic parameters statistics for cameras [' num2str(camInds) ']'];

subplot(2,5,1)
if size(RadialDistortion,2)==3
    boxplot(RadialDistortion,'Labels',{'k1','k2','k3'});
else
    boxplot(RadialDistortion,'Labels',{'k1','k2'});
end
ylim([-max(max(abs(RadialDistortion))) max(max(abs(RadialDistortion)))]);
hline = refline(0,0); hline.Color = 'g'; hline.LineStyle='--';
title('Radial Distorsion');
subplot(2,5,2,'align')
boxplot(TangentialDistortion,'Labels',{'p1','p2'});
%     ylim([-max(max(abs(TangentialDistortion))) max(max(abs(TangentialDistortion)))]);
hline =refline(0,0); hline.Color = 'g'; hline.LineStyle='--';
title('Tangential Distorsion');
subplot(2,5,3,'align')
boxplot(Skew,'Labels',{'s'});
hline = refline(0,0); hline.Color = 'g'; hline.LineStyle='--';
title('Skew');
subplot(2,5,4,'align')
boxplot(FocalLength,'Labels',{'fx','fy'});
title('Focal Length [pixel]');
subplot(2,5,5,'align')
boxplot(PrincipalPoint,'Labels',{'px','py'});
% plot principal points refline only if all cameras has same size
if sum(sum(imageSizeAll-mean(imageSizeAll)))==0
    hline =refline(0,imageSizeAll(1,1)/2); hline.Color = 'g'; hline.LineStyle='--';
    hline =refline(0,imageSizeAll(1,2)/2); hline.Color = 'g'; hline.LineStyle='--';
end
title('Principal Point [pixel]');

subplot(2,5,6,'align')
if size(RadialDistortionError,2)==3
    boxplot(RadialDistortionError,'Labels',{'k1','k2','k3'});
else
    boxplot(RadialDistortionError,'Labels',{'k1','k2'});
end
title('Radial Distorsion Error');
subplot(2,5,7,'align')
boxplot(TangentialDistortionError,'Labels',{'p1','p2'});
title('Tangential Distorsion Error');
subplot(2,5,8,'align')
boxplot(SkewError,'Labels',{'s'});
title('Skew Error');
subplot(2,5,9,'align')
boxplot(FocalLengthError,'Labels',{'fx','fy'});
title('Focal Length Error [pixel]');
subplot(2,5,10,'align')
boxplot(PrincipalPointError,'Labels',{'px','py'});
title('Principal Point Error [pixel]');

%% AFTER CORRECTION

tab(2) = uitab(tabgp,'Title','After distortion correction');
axes('Parent',tab(2));

for ic=1:Ncam
    RadialDistortion(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.RadialDistortion;
    RadialDistortionError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.RadialDistortionError;
    TangentialDistortion(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.TangentialDistortion;
    TangentialDistortionError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.TangentialDistortionError;
    Skew(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.Skew;
    SkewError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.SkewError;
    FocalLength(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.FocalLength;
    FocalLengthError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.FocalLengthError;
    PrincipalPoint(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.PrincipalPoint;
    PrincipalPointError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.PrincipalPointError;
    imageSizeAll(ic,:)=cameraCBparametersAllCams{ic}.imagesInfo.imageSize;
    camInds(ic)=cameraCBparametersAllCams{ic}.icam;
end

subplot(2,5,1)
if size(RadialDistortion,2)==3
    boxplot(RadialDistortion,'Labels',{'k1','k2','k3'});
else
    boxplot(RadialDistortion,'Labels',{'k1','k2'});
end
ylim([-max(max(abs(RadialDistortion))) max(max(abs(RadialDistortion)))]);
hline = refline(0,0); hline.Color = 'g'; hline.LineStyle='--';
title('Radial Distorsion');
subplot(2,5,2,'align')
boxplot(TangentialDistortion,'Labels',{'p1','p2'});
%     ylim([-max(max(abs(TangentialDistortion))) max(max(abs(TangentialDistortion)))]);
hline =refline(0,0); hline.Color = 'g'; hline.LineStyle='--';
title('Tangential Distorsion');
subplot(2,5,3,'align')
boxplot(Skew,'Labels',{'s'});
hline = refline(0,0); hline.Color = 'g'; hline.LineStyle='--';
title('Skew [deg]');
subplot(2,5,4,'align')
boxplot(FocalLength,'Labels',{'fx','fy'});
title('Focal Length [pixel]');
subplot(2,5,5,'align')
boxplot(PrincipalPoint,'Labels',{'px','py'});
% plot principal points refline only if all cameras has same size
if sum(sum(imageSizeAll-mean(imageSizeAll)))==0
    hline =refline(0,imageSizeAll(1,1)/2); hline.Color = 'g'; hline.LineStyle='--';
    hline =refline(0,imageSizeAll(1,2)/2); hline.Color = 'g'; hline.LineStyle='--';
end
title('Principal Point [pixel]');

subplot(2,5,6,'align')
if size(RadialDistortionError,2)==3
    boxplot(RadialDistortionError,'Labels',{'k1','k2','k3'});
else
    boxplot(RadialDistortionError,'Labels',{'k1','k2'});
end
title('Radial Distorsion Error');
subplot(2,5,7,'align')
boxplot(TangentialDistortionError,'Labels',{'p1','p2'});
title('Tangential Distorsion Error');
subplot(2,5,8,'align')
boxplot(SkewError,'Labels',{'s'});
title('Skew Error [deg]');
subplot(2,5,9,'align')
boxplot(FocalLengthError,'Labels',{'fx','fy'});
title('Focal Length Error [pixel]');
subplot(2,5,10,'align')
boxplot(PrincipalPointError,'Labels',{'px','py'});
title('Principal Point Error [pixel]');

%% BEFORE AND AFTYER

tab(3) = uitab(tabgp,'Title','Before (b) and after (a) distortion correction');
axes('Parent',tab(3));

for ic=1:Ncam
    RadialDistortion(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.RadialDistortion;
    RadialDistortionError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.RadialDistortionError;
    TangentialDistortion(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.TangentialDistortion;
    TangentialDistortionError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.TangentialDistortionError;
    Skew(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.Skew;
    SkewError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.SkewError;
    FocalLength(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.FocalLength;
    FocalLengthError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.FocalLengthError;
    PrincipalPoint(ic,:)=cameraCBparametersAllCams{ic}.cameraParameters.PrincipalPoint;
    PrincipalPointError(ic,:)=cameraCBparametersAllCams{ic}.estimationErrors.IntrinsicsErrors.PrincipalPointError;
    imageSizeAll(ic,:)=cameraCBparametersAllCams{ic}.imagesInfo.imageSize;

    RadialDistortionAUD(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.RadialDistortion;
    RadialDistortionErrorAUD(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.RadialDistortionError;
    TangentialDistortionAUD(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.TangentialDistortion;
    TangentialDistortionErrorAUD(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.TangentialDistortionError;
    SkewAUD(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.Skew;
    SkewErrorAUD(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.SkewError;
    FocalLengthAUD(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.FocalLength;
    FocalLengthErrorAUD(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.FocalLengthError;
    PrincipalPointAUD(ic,:)=cameraCBparametersAllCams{ic}.cameraParametersAUD.PrincipalPoint;
    PrincipalPointErrorAUD(ic,:)=cameraCBparametersAllCams{ic}.estimationErrorsAUD.IntrinsicsErrors.PrincipalPointError;

    camInds(ic)=cameraCBparametersAllCams{ic}.icam;
end

subplot(2,14,[1 2 3])
if size(RadialDistortion,2)==3
    boxplot([RadialDistortion(:,1) RadialDistortionAUD(:,1) RadialDistortion(:,2) RadialDistortionAUD(:,2) RadialDistortion(:,3) RadialDistortionAUD(:,3)],'Labels',{'k1(b)','k1(a)','k2(b)','k2(a)','k3(b)','k3(a)'});
    ylim([-max(max(abs([RadialDistortion RadialDistortionAUD]))) max(max(abs([RadialDistortion RadialDistortionAUD])))]);
    ylimits=ylim;
    line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');
    line([4.5 4.5],[ylimits(1) ylimits(2)],'color','k');
else
    boxplot([RadialDistortion(:,1) RadialDistortionAUD(:,1) RadialDistortion(:,2) RadialDistortionAUD(:,2)],'Labels',{'k1(b)','k1(a)','k2(b)','k2(a)'});
    ylim([-max(max(abs([RadialDistortion RadialDistortionAUD]))) max(max(abs([RadialDistortion RadialDistortionAUD])))]);
    ylimits=ylim;
    line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');
end
hline = refline(0,0); hline.Color = 'g'; hline.LineStyle='--';
title('Radial Distorsion Parameters');

subplot(2,14,[5 6],'align')
boxplot([TangentialDistortion(:,1) TangentialDistortionAUD(:,1) TangentialDistortion(:,2) TangentialDistortionAUD(:,2)],'Labels',{'p1(b)','p1(a)','p2(b)','p2(a)'});
ylimits=ylim;
line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');
%     ylim([-max(max(abs(TangentialDistortion))) max(max(abs(TangentialDistortion)))]);
hline =refline(0,0); hline.Color = 'g'; hline.LineStyle='--';
title('Tangential Distorsion Parameters');

subplot(2,14,8,'align')
boxplot([Skew SkewAUD],'Labels',{'s(b)','s(a)'});
hline = refline(0,0); hline.Color = 'g'; hline.LineStyle='--';
title('Skew Parameter');

subplot(2,14,[10 11],'align')
boxplot([FocalLength(:,1) FocalLengthAUD(:,1) FocalLength(:,2) FocalLengthAUD(:,2)],'Labels',{'fx(b)','fx(a)','fy(b)','fy(a)'});
ylimits=ylim;
line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');
title('Focal Length [pixel]');

subplot(2,14,[13 14],'align')
boxplot([PrincipalPoint(:,1) PrincipalPointAUD(:,1) PrincipalPoint(:,2) PrincipalPointAUD(:,2)],'Labels',{'px(b)','px(a)','py(b)','py(a)'});
ylimits=ylim;
line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');

% plot principal points refline only if all cameras has same size
if sum(sum(imageSizeAll-mean(imageSizeAll)))==0
    hline =refline(0,imageSizeAll(1,1)/2); hline.Color = 'g'; hline.LineStyle='--';
    hline =refline(0,imageSizeAll(1,2)/2); hline.Color = 'g'; hline.LineStyle='--';
end
title('Principal Point [pixel]');

subplot(2,14,[15 16 17],'align')
if size(RadialDistortionError,2)==3
    boxplot([RadialDistortionError(:,1) RadialDistortionErrorAUD(:,1) RadialDistortionError(:,2) RadialDistortionErrorAUD(:,2) RadialDistortionError(:,3) RadialDistortionErrorAUD(:,3)],'Labels',{'k1(b)','k1(a)','k2(b)','k2(a)','k3(b)','k3(a)'})
    ylimits=ylim;
    line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');
    line([4.5 4.5],[ylimits(1) ylimits(2)],'color','k');
else
    boxplot([RadialDistortionError(:,1) RadialDistortionErrorAUD(:,1) RadialDistortionError(:,2) RadialDistortionErrorAUD(:,2)],'Labels',{'k1(b)','k1(a)','k2(b)','k2(a)'})
    ylimits=ylim;
    line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');
end
title('Radial Distorsion Errors');

subplot(2,14,[19 20],'align')
boxplot([TangentialDistortionError(:,1) TangentialDistortionErrorAUD(:,1) TangentialDistortionError(:,2) TangentialDistortionErrorAUD(:,2)],'Labels',{'p1(b)','p1(a)','p2(b)','p2(a)'});
ylimits=ylim;
line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');
title('Tangential Distorsion Errors');

subplot(2,14,22,'align')
boxplot([SkewError SkewErrorAUD],'Labels',{'s(b)','s(a)'});
title('Skew Error');

subplot(2,14,[24 25],'align')
boxplot([FocalLengthError(:,1) FocalLengthErrorAUD(:,1) FocalLengthError(:,2) FocalLengthErrorAUD(:,2)],'Labels',{'fx(b)','fx(a)','fy(b)','fy(a)'});
ylimits=ylim;
line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');
title('Focal Length Errors [pixel]');

subplot(2,14,[27 28],'align')
boxplot([PrincipalPointError(:,1) PrincipalPointErrorAUD(:,1) PrincipalPointError(:,2) PrincipalPointErrorAUD(:,2)],'Labels',{'px(b)','px(a)','py(b)','py(a)'});
ylimits=ylim;
line([2.5 2.5],[ylimits(1) ylimits(2)],'color','k');
title('Principal Point Errors [pixel]');

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
%
% generate training data
function [dres_train, dres_det, labels] = generate_training_data(seq_idx, dres_image, opt)

% is_show = 0;

seq_name = opt.mot2d_train_seqs{seq_idx};
seq_set = 'train';

% read detections
filename = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'det', 'det.txt');
dres_det = read_mot2dres(filename);

% read ground truth
filename = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'gt', 'gt.txt');
dres_gt = read_mot2dres(filename);
dres_gt = fix_groundtruth(seq_name, dres_gt);
y_gt = dres_gt.y + dres_gt.h;

% collect true positives and false alarms from detections
num = numel(dres_det.fr);
labels = zeros(num, 1);
overlaps = zeros(num, 1);
for i = 1:num
    fr = dres_det.fr(i);
    index = find(dres_gt.fr == fr);
    if isempty(index) == 0
        overlap = calc_overlap(dres_det, i, dres_gt, index);
        o = max(overlap);
        if o < opt.overlap_neg
            labels(i) = -1;
        elseif o > opt.overlap_pos
            labels(i) = 1;
        else
            labels(i) = 0;
        end
        overlaps(i) = o;
    else
        overlaps(i) = 0;
        labels(i) = -1;
    end
end

% build the training sequences
ids = unique(dres_gt.id);
dres_train = [];
count = 0;
for i = 1:numel(ids)
    index = find(dres_gt.id == ids(i));
    dres = sub(dres_gt, index);

    % check if the target is occluded or not
    num = numel(dres.fr);
    dres.occluded = zeros(num, 1);
    dres.covered = zeros(num, 1);
    dres.overlap = zeros(num, 1);
    dres.r = zeros(num, 1);
    dres.area_inside = zeros(num, 1);
    y = dres.y + dres.h;
    for j = 1:num
        fr = dres.fr(j);
        index = find(dres_gt.fr == fr & dres_gt.id ~= ids(i));

        if isempty(index) == 0
            [~, ov] = calc_overlap(dres, j, dres_gt, index);
            ov(y(j) > y_gt(index)) = 0;
            dres.covered(j) = max(ov);
        end

        if dres.covered(j) > opt.overlap_occ
            dres.occluded(j) = 1;
        end

        % overlap with detections
        index = find(dres_det.fr == fr);
        if isempty(index) == 0
            overlap = calc_overlap(dres, j, dres_det, index);
            [o, ind] = max(overlap);
            dres.overlap(j) = o;
            dres.r(j) = dres_det.r(index(ind));

            % area inside image
            [~, overlap] = calc_overlap(dres_det, index(ind), dres_image, fr);
            dres.area_inside(j) = overlap;
        end
    end

    % start with bounding overlap > opt.overlap_pos and non-occluded box
    index = find(dres.overlap > opt.overlap_pos & dres.covered == 0 & dres.area_inside > opt.exit_threshold);
    if isempty(index) == 0
        index_start = index(1);
        count = count + 1;
        dres_train{count} = sub(dres, index_start:num);
    end

    % show gt
%      if is_show
%         for j = 1:numel(dres_train{i}.fr)
%             fr = dres_train{i}.fr(j);
%             filename = fullfile(opt.mot, opt.mot2d, seq_set, seq_name, 'img1', sprintf('%06d.jpg', fr));
%             disp(filename);
%             I = imread(filename);
%             figure(1);
%             show_dres(fr, I, 'GT', dres_train{i});
%             pause;
%         end
%     end
end

% handle occlusion by a pole in PETS09-S2L1
if strcmp(seq_name, 'PETS09-S2L1') == 1
    dres_train([9, 15, 19]) = [];
    pole = [409.0000, 102.0000, 39.0000, 294.0000];
    dres_pole.x = pole(1);
    dres_pole.y = pole(2);
    dres_pole.w = pole(3);
    dres_pole.h = pole(4);
    for i = 1:numel(dres_train)
        dres = dres_train{i};
        [~, ~, overlap] = calc_overlap(dres_pole, 1, dres, 1:numel(dres.fr));
        dres.covered = max([dres.covered, overlap'], [], 2);
        dres.occluded(dres.covered > opt.overlap_occ) = 1;
        dres_train{i} = dres;
    end
end

fprintf('%s: %d positive sequences\n', seq_name, numel(dres_train));

author: yuxng


text:
function z = rel_entr( x, y )

%REL_ENTR   Scalar relative entropy.
%   REL_ENTR(X,Y) returns an array of the same size as X+Y with the
%   relative entropy function applied to each element:
%                      { X.*LOG(X./Y) if X >  0 & Y >  0,
%      REL_ENTR(X,Y) = { 0            if X == 0 & Y >= 0,
%                      { +Inf         otherwise.
%   X and Y must either be the same size, or one must be a scalar. If X and
%   Y are vectors, then SUM(REL_ENTR(X,Y)) returns their relative entropy.
%   If they are PDFs (that is, if X>=0, Y>=0, SUM(X)==1, SUM(Y)==1) then
%   this is equal to their Kullback-Liebler divergence SUM(KL_DIV(X,Y)).
%   -SUM(REL_ENTR(X,1)) returns the entropy of X.
%
%   Disciplined convex programming information:
%       REL_ENTR(X,Y) is convex in both X and Y, nonmonotonic in X, and
%       nonincreasing in Y. Thus when used in CVX expressions, X must be
%       real and affine and Y must be concave. The use of REL_ENTR(X,Y) in
%       an objective or constraint will effectively constrain both X and Y
%       to be nonnegative, hence there is no need to add additional
%       constraints X >= 0 or Y >= 0 to enforce this.

error(nargchk(2,2,nargin));
if ~isreal( x ) || ~isreal( y ),
    error( 'Arguments must be real.' );
end
t1 = x < 0  | y <= 0;
t2 = x == 0 & y >= 0;
x  = max( x, realmin );
y  = max( y, realmin );
z  = x .* log( x ./ y );
z( t1 ) = +Inf;
z( t2 ) = 0;

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
% normalize a vector
function result = norma(v)
    result = v / norm(v);
end

author: layumi

repo: Image-Text-Embedding


text:
 function y = mtimes(ob, x)
%function y = mtimes(ob, x)
% y = G * x	or x = G' * y
% Copyright 2002-2-20, Jeff Fessler, The University of Michigan

%
% scalar * G
%
if isa(ob, 'double') & length(x) == 1 & isa(x, 'Gtomo2_dsc')
	y = x;
	y.scale = ob;
	return
end

if ob.apower ~= 1, error notdone, end

%
% full forward projection
%
if ~ob.is_transpose

	% if needed, expand concise column(s)
	if ob.is_masked
		idim = size(x);
		np = sum(ob.mask(:));
		nxy = numel(ob.mask);
		if idim(1) ~= np
			error 'size mismatch'
		end
		x = embed(x, ob.mask);
		x = reshape(x, nxy, idim(2));
	end

	y = wtfmex('dsc,proj', ob.arg', single(x), uint8(ob.mask), ...
		int32(0), int32(1), int32(ob.nthread), int32(ob.chat));

%
% full back-projection
%
else
	y = wtfmex('dsc,back', ob.arg', single(x), uint8(ob.mask), ...
		int32(0), int32(1), int32(ob.nthread), int32(ob.chat));

	if ob.is_masked
		if size(y,1) == numel(ob.mask)
			y = y(ob.mask,:);	% [nxy,nz] -> [np,nz]
		else
			y = y(ob.mask);		% [nx,ny] -> [np]
		end
	end
end

y = ob.scale * double(y);	% trick: no '(:)' here; let wtfmex do the job!

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function result = cached_argument(name,default)
% Helper: return an arg_nogui() specifier with the given name/default value
% Result = cached_argument(Name,Default)
%
% This function is an optimization for internal use by the argument system.
%
% In:
%   Name : name of the argument
%
%   Default : default value of the argument
%
% Out:
%   Result : an argument specifier struct, with value pre-assigned.
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2014-02-03

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

persistent keys values;
try
    key = [name '_' char('a'+default)];
    try
        result = values{strcmp(keys,key)};
    catch %#ok<CTCH>
        result = arg_nogui(name,default);
        result = feval(result{1},[],result{2}{:});
        result.value = result.defaults{1};
        if ~iscell(keys)
            keys = {key};
            values = {result};
        else
            keys{end+1} = key;
            values{end+1} = result;
        end
    end
catch %#ok<CTCH>
    result = arg_nogui(name,default);
    result = feval(result{1},[],result{2}{:});
    result.value = result.defaults{1};
end


text:
classdef CatalogManager < handle
    % Manage a bunch of catalogs, filtering according to each catalog's associated filter function
    %
    % Usage:
    %   % associate the CatalogManager with a catalog
    %   cm = CatalogManger(someZmapCatalog)
    %
    %   % add a 'live' subsetter that will track the large (M>5) events
    %
    %   cm.AddSubset( 'big_events', @(c)c.Magnitude > 5)
    %
    %   % add another subsetter that will track shallow northen hemisphere events
    %
    %   cm.AddSubset( 'shallow_northern', @(c)c.Depth<5 & c.Latitude>0)
    %
    %
    %   % if you would like ot modify the filter, then use 'ChangeFilter'
    %   %     - using AddSubset would be an error, since it already exists
    %   cm.ChangeFilter(('shallow_northern', @(c)c.Depth<8 & c.Latitude>0)
    %
    %
    %   % to get the subcatalog, use 'GetSubset'
    %   shallow_northern_cat = cm.GetSubset('shallow_northern')
    %
    %   % change the catalog and get some new results...
    %   cm.RawCatalog = someOtherZmapCatalog
    %   shallow_northern_cat = cm.GetSubset('shallow_northern')
    %

    properties
        RawCatalog % this is the catalog that will be filtered
    end

    properties(Access=private)
        Catalogs        containers.Map
    end

    events
        RawCatalogChanged
    end

    methods
        function obj = CatalogManager(catalog)
            obj.RawCatalog = catalog;
        end

        function AddSubset(obj, Tag, filterFcn, ~)
            subcatalog.FilterFcn = filterFcn; % filterFcn returns a logical mask
            subcatalog.Filter = subcatalog.FilterFcn(obj.RawCatalog);
            if ismember(Tag, keys(obj.Catalogs))
                error('%s is already being used. Did you mean to CHANGE the filter instead?',Tag)
            end
            obj.Catalogs(Tag) = subcatalog;
        end

        function catalog = GetSubset(obj, Tag)
            subcatalog = obj.Catalogs(Tag);
            if numel(subcatalog.Filter) ~= obj.RawCatalog.Count
                subcatalog.Filter = subcatalog.FilterFcn(obj.RawCatalog);
            end
            catalog = obj.RawCatalog.subset(subcatalog.Filter);
        end

        function filt = GetFilter(obj, Tag, position)
            % position can be 'first' or 'last' or empty (empty returns all)
            subcatalog = obj.Catalogs(Tag);
            if numel(subcatalog.Filter) ~= obj.RawCatalog.Count
                subcatalog.Filter = subcatalog.FilterFcn(obj.RawCatalog);
            end
            filt = subcatalog.Filter;
            if exist('position','var')
            	p = find(filt, 1, position);
                filt = false(size(filt));
                filt(p) = true;
            end
        end

        function ChangeFilter(obj, Tag, filterFcn)
            subcatalog.FilterFcn = filterFcn; % filterFcn returns a logical mask
            subcatalog.Filter = subcatalog.FilterFcn(obj.RawCatalog);
            obj.Catalogs(Tag) = subcatalog;
        end

        function RemoveSubset(obj, Tag)
            remove(obj.Catalogs,Tag);
        end

        function RecalculateFilter(obj, Tag)
            subcatalog = obj.Catalogs(Tag);
            subcatalog.Filter = subcatalog.FilterFcn(obj.RawCatalog);
            obj.Catalogs(Tag) = subcatalog;
        end

        function recalculateAllFilters(obj)
            for key = keys(obj.Catalogs)
                obj.RecalculateFilter(key{1});
            end
        end

        function disp(obj)
            getIndentedSummary = @(ca,level) "     "+ strrep(summary(ca,level),newline,[newline, '       ']);
            disp('CatalogManger with:')
            disp('  RawCatalog:')
            disp(getIndentedSummary(obj.RawCatalog,'simple'));
            fprintf('\n\nWith Tagged entries:\n')
            for key = keys(obj.Catalogs)
                Tag = key{1};
                c = obj.GetSubset(Tag);
                entry = obj.Catalogs(Tag);
                fprintf("  %-10s : %s  yielding\n","'"+Tag+"'", func2str(entry.FilterFcn));
                disp(getIndentedSummary(c,'simple'));
                fprintf('\n\n');
            end
        end
    end
end

author: CelsoReyes


text:
function [xs,nxs] = get_position_and_orientation_ls(xs,conf);
%GET_POSITION_AND_ORIENTATION_LS position and orientation for a line source
%
%   Usage: [xs,nxs] = get_position_and_orientation_ls(xs,conf);
%
%   Input parameters:
%       xs          - combined position and orientation / m [nx3] or [nx6]
%       args        - list of args
%
%   Output parameters:
%       xs          - position of line source / m [nx3]
%       nxs         - orientation of line source / m [nx3]
%
%   GET_POSITION_AND_ORIENTATION_LS(xs,conf) returns the position and
%   orientation of a point source. The orientation is a vector perpendicular to
%   the traveling direction of the line source. For 2D or 2.5D the orientation
%   will always be returned as [0 0 1]. For 3D [0 0 1] will be returned if no
%   explicit orientation is given.
%
%   See also: secondary_source_selection, driving_function_mono_wfs

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 2;
nargmax = 2;
narginchk(nargmin,nargmax);
isargnumeric(xs);
isargstruct(conf);

%% ===== Configuration ===================================================
dimension = conf.dimension;

%% ===== Checking for vector =============================================
% Handling of line source orientation
if (strcmp('2D',dimension) ||  strcmp('2.5D',dimension))
    % Ignore orientation for 2D and 2.5D
    if size(xs,2)>3
        warning('%s: %s-WFS ignores virtual line source orientation.', ...
            upper(mfilename),dimension);
    end
    xs(:,3) = 0;
    nxs = repmat([0 0 1],[size(xs,1) 1]);
else
    if size(xs,2)~=6
        warning('%s: set ls orientation to [0 0 1]',upper(mfilename));
        nxs = repmat([0 0 1],[size(xs,1) 1]);
    else
        nxs = xs(:,4:6) / norm(xs(1,4:6),2);
    end
end
xs = xs(:,1:3);

author: sfstoolbox

repo: sfs-matlab


text:
function res = chowlin_co_W(Y,x,ta,sc,opC)
% PURPOSE: Temporal disaggregation using the Chow-Lin method
%         (quarterly rho derived from Cochrane-Orcutt annual rho)
%          Without pretesting for intercept
% ------------------------------------------------------------
% SYNTAX: res = chowlin_co_W(Y,x,ta,sc,opC);
% ------------------------------------------------------------
% OUTPUT: res: a structure
%           res.meth    ='Chow-Lin';
%           res.ta      = type of disaggregation
%           res.type    = method of estimation
%           res.N       = nobs. of low frequency data
%           res.n       = nobs. of high-frequency data
%           res.pred    = number of extrapolations
%           res.sc       = frequency conversion between low and high freq.
%           res.p       = number of regressors (including intercept)
%           res.Y       = low frequency data
%           res.x       = high frequency indicators
%           res.y       = high frequency estimate
%           res.y_dt    = high frequency estimate: standard deviation
%           res.y_lo    = high frequency estimate: sd - sigma
%           res.y_up    = high frequency estimate: sd + sigma
%           res.u       = high frequency residuals
%           res.U       = low frequency residuals
%           res.beta    = estimated model parameters
%           res.beta_sd = estimated model parameters: standard deviation
%           res.beta_t  = estimated model parameters: t ratios
%           res.rho     = innovational parameter
%           res.aic     = Information criterion: AIC
%           res.bic     = Information criterion: BIC
%           res.co      = Cochrane-Orcutt regression (see LeSage
%           Econometric Toolbox)
% ------------------------------------------------------------
% INPUT: Y: Nx1 ---> vector of low frequency data
%        x: nxp ---> matrix of high frequency indicators (without intercept)
%        ta: type of disaggregation
%            ta=1 ---> sum (flow)
%            ta=2 ---> average (index)
%        sc: number of high frequency data points for each low frequency data points
%            sc= 3 ---> quarterly to monthly
%            sc= 4 ---> annual to quarterly
%        opC: 1x1 option related to intercept
%            opc = 0 : no intercept in hf model
%            opc = 1 : intercept in hf model
% ------------------------------------------------------------
% LIBRARY: aggreg, olsc
% ------------------------------------------------------------
% SEE ALSO: chowlin, litterman, fernandez, td_plot, td_print
% ------------------------------------------------------------
% REFERENCE: Chow, G. and Lin, A.L. (1971) "Best linear unbiased
% distribution and extrapolation of economic time series by related
% series", Review of Economic and Statistics, vol. 53, n. 4, p. 372-375.
% Bournay, J. y Laroque, G. (1979) "Reflexions sur la methode d'elaboration
% des comptes trimestriels", Annales de l'INSEE, n. 36, p. 3-30.

% written by:
%  Enrique M. Quilis
%  Macroeconomic Research Department
%  Ministry of Economy and Competitiveness
%  <enrique.quilis@mineco.es>

% Version 1.1 [August 2006]

t0=clock;

% ------------------------------------------------------------
% Checks

if ((sc ~= 3) & (sc ~= 4))
    error ('*** THE FREQUENCY CONVERSION SHOULD BE 3 OR 4 ***');
 end

if ((ta == 3) | (ta == 4))
    error ('*** THIS RELEASE DOES NOT PERFORM INTERPOLATION ***');
end

% ------------------------------------------------------------
% Size of the problem

[N,M] = size(Y);    % Size of low-frequency input
[n,p] = size(x);    % Size of p high-frequency inputs (without intercept)

% ------------------------------------------------------------
% Preparing the X matrix: including an intercept if opC==1

if (opC == 1)
   e=ones(n,1);
   x=[e x];       % Expanding the regressor matrix
   p=p+1;         % Number of p high-frequency inputs (plus intercept)
end

% ------------------------------------------------------------
% Generating the aggregation matrix

C = aggreg(ta,N,sc);

% -----------------------------------------------------------
% Expanding the aggregation matrix to perform
% extrapolation if needed.

if (n > sc * N)
   pred=n-sc*N;           % Number of required extrapolations
   C=[C zeros(N,pred)];
else
   pred=0;
end

% -----------------------------------------------------------
% Temporal aggregation of the indicators

X=C*x;

% ------------------------------------------------------------
% Computing annual rho by means of Cochrane-Orcutt

rex = olsc(Y,X); %Econometric Toolbox (LeSage, 1999)

Ra=rex.rho; %Getting annual rho

if (Ra < 0)
    error ('*** ANNUAL RHO IS NEGATIVE. END OF PROGRAM ***');
end

if (Ra >= 1)
   Ra = 0.99;
end

r = 0:0.01:0.99; nr = length(r); R=ones(nr,1);
switch sc
case 3 %Quarterly to monthly
   for i=1:nr
      R(i) = (r(i)^5+2*r(i)^4+3*r(i)^3+2*r(i)^2+r(i)) / (2*r(i)^2+4*r(i)+3);
   end
case 4 %Annual to quarterly
   % ------------------------------------------------------------
   % Evaluating the function that relates annual and quarterly rho
   for i=1:nr
      R(i) = (r(i)*(r(i)+1)*(r(i)^2+1)^2) / (2*(r(i)^2+r(i)+2));
   end
end

% -----------------------------------------------------------
% Determination of optimal rho

[aux,h] = min(abs(Ra-R));
rho = r(h);

% -----------------------------------------------------------
% Final estimation with optimal rho

I=eye(n); w=I;
LL = diag(-ones(n-1,1),-1);

Aux=I+rho*LL;
Aux(1,1)=sqrt(1-rho^2);
w=inv(Aux'*Aux);           % High frequency VCV matrix (without sigma_a)
W=C*w*C';                  % Low frequency VCV matrix (without sigma_a)
Wi=inv(W);
beta=(X'*Wi*X)\(X'*Wi*Y);  % beta estimator
U=Y-X*beta;                % Low frequency residuals
scp=U'*Wi*U;               % Weighted least squares
sigma_a=scp/(N-p);         % sigma_a estimator
L=w*C'*Wi;                 % Filtering matrix
u=L*U;                     % High frequency residuals

% -----------------------------------------------------------
% Temporally disaggregated time series

y = x*beta + u;

% -----------------------------------------------------------
% Information criteria
% Note: p is expanded to include the innovational parameter

aic=log(sigma_a)+2*(p+1)/N;
bic=log(sigma_a)+log(N)*(p+1)/N;

% -----------------------------------------------------------
% VCV matrix of high frequency estimates

sigma_beta=sigma_a*inv(X'*Wi*X);

VCV_y=sigma_a*(eye(n)-L*C)*w+(x-L*X)*sigma_beta*(x-L*X)';

d_y=sqrt((diag(VCV_y)));   % Std. dev. of high frequency estimates
y_li=y-d_y;           % Lower lim. of high frequency estimates
y_ls=y+d_y;           % Upper lim. of high frequency estimates

% -----------------------------------------------------------
% -----------------------------------------------------------
% Loading the structure

res.meth='Chow-Lin';

% -----------------------------------------------------------
% Basic parameters

res.ta        = ta;
res.N         = N;
res.n         = n;
res.pred      = pred;
res.sc        = sc;
res.p         = p;
res.type      = 3;  % For output convenience
res.opC       = opC;

% -----------------------------------------------------------
% Series

res.Y         = Y;
res.x         = x;
res.y         = y;
res.y_dt      = d_y;
res.y_lo      = y_li;
res.y_up      = y_ls;

% -----------------------------------------------------------
% Residuals

res.u         = u;
res.U         = U;

% -----------------------------------------------------------
% Parameters

res.beta      = beta;
res.beta_sd   = sqrt(diag(sigma_beta));
res.beta_t    = beta./sqrt(diag(sigma_beta));
res.rho       = rho;

% -----------------------------------------------------------
% Information criteria

res.aic       = aic;
res.bic       = bic;

% -----------------------------------------------------------
% Cochrane-Orcutt results

res.co       = rex;

% -----------------------------------------------------------
% Objective function

res.val       = R;
res.r         = r;
res.valA      = Ra * ones(nr,1);
res.ropt      = rho * ones(nr,1);

% -----------------------------------------------------------
% Elapsed time

res.et        = etime(clock,t0);


text:
function mtrExportFibersFromMatrix(pathwayMatrix, pathwayLengthVec, statsVec, statsNames, fileName, xformToMm)
%Exports pathway matrix (3,max_path_length,num_paths) to .pdb format
%
% This routine does not appear to be used much.
%
%   mtrExportFibersFromMatrix(pathwayMatrix, pathwayLengthVec, fileName, xformToAcPc, xform)
%
% See also: mtrNotes
%
% Examples
%
%
% HISTORY:
% 2008.01.13 AJS: wrote it.
%
%
% Stanford VISTA Team

if(~exist('xformToMm','var')),  xformToMm = eye(4); end

% Create blank database for finding some of the sizes to write out
pdb = mtrPathwayDatabase();

% PDB structure will have stat header info, but fg will retain the actual
% statistics per path to record
mmPerStep = sqrt(sum((pathwayMatrix(:,1,1)-pathwayMatrix(:,2,1)).^2));
[pdb,fg] = makeStatsFromStatsVec(pdb,pathwayLengthVec,mmPerStep, statsVec, statsNames);
numstats = length( pdb.pathway_statistic_headers );
numalgs = 0;

% Redo offset
offset = 4*4 + 16*8 + numstats * getStatisticsHeaderSizeWordAligned(pdb) + numalgs * getAlgorithmsHeaderSizeWordAligned(pdb); % 4 uint, 16 doubles

fid = fopen(fileName, 'wb');

% Writing header
fwrite(fid,offset,'uint');
xformToMm = xformToMm';
fwrite(fid,xformToMm(:),'double');
fwrite(fid,numstats,'uint');

% Write stat header
saveStatisticsHeader(pdb,fid);

fwrite(fid,numalgs,'uint');
% Write algorithm header
%saveAlgorithmHeader(pdb,fid);

% Write version info
version = 1;
fwrite(fid,version,'uint');

%Writing path info
numpaths = length(pathwayLengthVec);
fwrite(fid,numpaths,'uint');
fileOffsets = zeros(1,numpaths);
for ff = 1:numpaths
    % Record current file position to write later
    fileOffsets(ff) = ftell(fid);

    path_offset = 3 * 4  + numstats*8; % 4 int, numstats double
    fwrite(fid,path_offset,'int');
    nn = pathwayLengthVec(ff);
    fwrite(fid,nn,'int');
    fwrite(fid,0,'int'); % algo type HACK
    fwrite(fid,1,'int'); % seedpointindex HACK ???

    % Write stats per path
    for as = 1:numstats
        fwrite(fid,fg.params{as}.stat(ff),'double');
    end

    % Writing path nodes
    %pos = fg.fibers{ff}.coords;
    pos = pathwayMatrix(:,1:pathwayLengthVec(ff),ff);
    %pos = mrAnatXformCoords(xformFromAcpc, pos')';
    % Change from 1 based to 0 based positions in mm space
    %pos = pos - repmat([2,2,2]',1,size(pos,2));
    fwrite(fid,pos,'double');

    % Writing stats values per position
    for as = 1:numstats
        % Just write garbage as fiber groups can't handle per point stats
        if( pdb.pathway_statistic_headers(as).is_computed_per_point )
            error('Cannot export fiber group with per point statistics!');
            %fwrite(fid,this.pathways(p).point_stat_array(as,:),'double');
        end
    end
end

% Database Footer
fwrite(fid,fileOffsets,'uint32');

fclose(fid);
return;

function [pdb, fg] = makeStatsFromStatsVec(pdb, lengthVec, mmPerStep, statsVec, statsNames)

fg = dtiNewFiberGroup();
fg.params = {};
fg.params{1}.name = 'Length';
for ff=1:length(lengthVec)
    fg.params{1}.stat(ff) = mmPerStep*(lengthVec(ff)-1);
end
s=length(pdb.pathway_statistic_headers)+1;
pdb = addStatisticHeader(pdb,'Length','NA', 1, 0, 1, 8877+s);

for ii=2:min(length(statsNames)+1,3)
    % Add up to two stats from the provided vectors
    fg.params{ii}.name = statsNames{ii-1};
    for ff=1:size(statsVec,2)
        fg.params{ii}.stat(ff) = statsVec(ii-1,ff);
    end
    s=length(pdb.pathway_statistic_headers)+1;
    pdb = addStatisticHeader(pdb,statsNames{ii-1},'NA', 1, 0, 1, 8877+s);
end

for ii=(length(statsNames)+2):3
    % Add blank stats to fill up to 3 required by pdb format.
    pdb = addBlankStat(pdb,'Blank');
    fg.params{ii}.name = 'Blank';
    fg.params{ii}.stat = zeros(1,length(lengthVec));
end

return;

function pdb = addBlankStat(pdb,agg_name)
s=length(pdb.pathway_statistic_headers)+1;
pdb = addStatisticHeader(pdb,agg_name,'NA', 1, 0, 1, 8877+s);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [logRS,logERS,V]=RSana(x,n,method,q)
%Syntax: [logRS,logERS,V]=RSana(x,n,method,q)
%____________________________________________
%
% Performs R/S analysis on a time series.
%
% logRS is the log(R/S).
% logERS is the Expectation of log(R/S).
% V is the V statistic.
% x is the time series.
% n is the vector with the sub-periods.
% method can take one of the following values
%  'Hurst' for the Hurst-Mandelbrot variation.
%  'Lo' for the Lo variation.
%  'MW' for the Moody-Wu variation.
%  'Parzen' for the Parzen variation.
% q can be either
%  a (non-negative) integer.
%  'auto' for the Lo's suggested value.
%
%
% References:
%
% Peters E (1991): Chaos and Order in the Capital Markets. Willey
%
% Peters E (1996): Fractal Market Analysis. Wiley
%
% Lo A (1991): Long term memory in stock market prices. Econometrica
% 59: 1279-1313
%
% Moody J, Wu L (1996): Improved estimates for Rescaled Range and Hurst
% exponents. Neural Networks in Financial Engineering, eds. Refenes A-P
% Abu-Mustafa Y, Moody J, Weigend A: 537-553, Word Scientific
%
% Hauser M (1997): Semiparametric and nonparametric testing for long
% memory: A Monte Carlo study. Empirical Economics 22: 247-271
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 1 Jan 2004.

if nargin<1 | isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
   % N is the time series length
   N=length(x);
end

if nargin<2 | isempty(n)==1
   n=1;
else
   % n must be either a scalar or a vector
   if min(size(n))>1
      error('n must be either a scalar or a vector.');
   end
   % n must be integer
   if n-round(n)~=0
       error('n must be integer.');
   end
   % n must be positive
   if n<=0
      error('n must be positive.');
   end
end

if nargin<4 | isempty(q)==1
   q=0;
else
    if q=='auto'
        t=autocorr(x,1);
        t=t(2);
        q=((3*N/2)^(1/3))*(2*t/(1-t^2))^(2/3);
    else
        % q must be a scalar
        if sum(size(q))>2
            error('q must be scalar.');
        end
        % q must be integer
        if q-round(q)~=0
            error('q must be integer.');
        end
        % q must be positive
        if q<0
            error('q must be positive.');
        end
    end
end

for i=1:length(n)

    % Calculate the sub-periods
    a=floor(N/n(i));

    % Make the sub-periods matrix
    X=reshape(x(1:a*n(i)),n(i),a);

    % Estimate the mean of each sub-period
    ave=mean(X);

    % Remove the mean from each sub-period
    cumdev=X-ones(n(i),1)*ave;

    % Estimate the cumulative deviation from the mean
    cumdev=cumsum(cumdev);

    % Estimate the standard deviation
    switch method
    case 'Hurst'
        % Hurst-Mandelbrot variation
        stdev=std(X);
    case 'Lo'
        % Lo variation
        for j=1:a
            sq=0;
            for k=0:q
                v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
                if k>0
                    sq=sq+(1-k/(q+1))*v(k+1);
                end
            end
            stdev(j)=sqrt(v(1)+2*sq);
        end
    case 'MW'
        % Moody-Wu variation
        for j=1:a
            sq1=0;
            sq2=0;
            for k=0:q
                v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
                if k>0
                    sq1=sq1+(1-k/(q+1))*(n(i)-k)/n(i)/n(i);
                    sq2=sq2+(1-k/(q+1))*v(k+1);
                end
            end
            stdev(j)=sqrt((1+2*sq1)*v(1)+2*sq2);
        end
    case 'Parzen'
        % Parzen variation
        if mod(q,2)~=0
            error('For the "Parzen" variation q must be dived by 2.');
        end
        for j=1:a
            sq1=0;
            sq2=0;
            for k=0:q
                v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
                if k>0 & k<=q/2
                    sq1=sq1+(1-6*(k/q)^2+6*(k/q)^3)*v(k+1);
                elseif k>0 & k>q/2
                    sq2=sq2+(1-(k/q)^3)*v(k+1);
                end
            end
            stdev(j)=sqrt(v(1)+2*sq1+2*sq2);
        end
    otherwise
        error('You should provide another value for "method".');
    end

    % Estiamte the rescaled range
    rs=(max(cumdev)-min(cumdev))./stdev;

    clear stdev

    % Take the logarithm of the mean R/S
    logRS(i,1)=log10(mean(rs));

    if nargout>1

        % Initial calculations fro the log(E(R/S))
        j=1:n(i)-1;
        s=sqrt((n(i)-j)./j);
        s=sum(s);

        % The estimation of log(E(R/S))
        logERS(i,1)=log10(s/sqrt(n(i)*pi/2));

        % Other estimations of log(E(R/S))
        %logERS(i,1)=log10((n(i)-0.5)/n(i)*s/sqrt(n(i)*pi/2));
        %logERS(i,1)=log10(sqrt(n(i)*pi/2));

    end

    if nargout>2
        % Estimate V
        V(i,1)=mean(rs)/sqrt(n(i));
    end

end


text:
function [cfg] = ft_spike_plot_isi(cfg, isih)

% FT_SPIKE_PLOT_ISI makes an inter-spike-interval bar plot.
%
% Use as
%   ft_spike_plot_isi(cfg, isih)
%
% Inputs:
%   ISIH is the output from FT_SPIKE_ISIHIST
%
% Configurations:
%   cfg.spikechannel     = string or index or logical array to to select 1 spike channel.
%                          (default = 1).
%   cfg.ylim             = [min max] or 'auto' (default)
%                          If 'auto', we plot from 0 to 110% of maximum plotted value);
%   cfg.plotfit          = 'yes' (default) or 'no'. This requires that when calling
%                          FT_SPIKESTATION_ISI, cfg.gammafit = 'yes'.
%
% Outputs:
%   hdl.fit              = handle for line fit. Use SET and GET to access.
%   hdl.isih             = handle for bar isi histogram. Use SET and GET to access.

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble provenance isih

% get the default options
cfg.spikechannel = ft_getopt(cfg,'spikechannel', 'all');
cfg.ylim         = ft_getopt(cfg,'ylim', 'auto');
cfg.plotfit      = ft_getopt(cfg,'plotfit', 'no');

% ensure that the options are valid
cfg = ft_checkopt(cfg,'spikechannel',{'cell', 'char', 'double'});
cfg = ft_checkopt(cfg,'ylim', {'char','ascendingdoublebivector'});
cfg = ft_checkopt(cfg,'plotfit', 'char', {'yes', 'no'});

% get the spikechannels
cfg.channel = ft_channelselection(cfg.spikechannel, isih.label);
spikesel    = match_str(isih.label, cfg.channel);
nUnits      = length(spikesel); % number of spike channels
if nUnits~=1, error('One spikechannel should be selected by means of cfg.spikechannel'); end

% plot the average isih
isiHdl = bar(isih.time,isih.avg(spikesel,:),'k');
set(isiHdl,'BarWidth', 1)

if strcmp(cfg.plotfit,'yes')

  if ~isfield(isih,'gamfit'), error('isih.gamfit should be present when cfg.plotfit = yes'); end

  % generate the probabilities according to the gamma model
  pGamma = gampdf(isih.time, isih.gamfit(spikesel,1), isih.gamfit(spikesel,2));
  sel    = isfinite(pGamma);
  pGamma = pGamma(sel)./sum(pGamma(sel));

  % scale the fit in the same way (total area is equal)
  sumIsih = sum(isih.avg(spikesel,:),2);
  pGamma  = pGamma*sumIsih;

  % plot the fit
  hold on
  fitHdl = plot(isih.time(sel), pGamma,'r');
else
  pGamma = 0;
end

try
  if strcmp(isih.cfg.outputunit, 'proportion')
    ylabel('probability')
  elseif strcmp(isih.cfg.outputunit, 'spikecount')
    ylabel('spikecount')
  end
catch
  ylabel('intensity');
end

xlabel('bin edges (sec)')

% set the x axis
set(gca,'XLim', [min(isih.time) max(isih.time)])

% set the y axis
if strcmp(cfg.ylim, 'auto')
  y = [isih.avg(spikesel,:) pGamma];
  cfg.ylim = [0 max(y(:))*1.1+eps];
end
set(gca,'YLim', cfg.ylim)
set(gca,'TickDir','out')

% store the handles
if strcmp(cfg.plotfit,'yes'),     H.fit    = fitHdl;    end
H.isih   = isiHdl;

% as usual, make sure that panning and zooming does not distort the y limits
set(zoom,'ActionPostCallback',{@mypostcallback,cfg.ylim,[min(isih.time) max(isih.time)]});
set(pan,'ActionPostCallback',{@mypostcallback,cfg.ylim,[min(isih.time) max(isih.time)]});

% do the general cleanup and bookkeeping at the end of the function

ft_postamble previous isih
ft_postamble provenance

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [] = mypostcallback(fig,evd,limY,limX)

% get the x limits and reset them
ax = evd.Axes;
xlim = get(ax(1), 'XLim');
ylim = get(ax(1), 'YLim');
if limX(1)>xlim(1), xlim(1) = limX(1); end
if limX(2)<xlim(2), xlim(2) = limX(2); end
if limY(1)>ylim(1), ylim(1) = limY(1); end
if limY(2)<ylim(2), ylim(2) = limY(2); end

set(ax(1), 'XLim',xlim)
set(ax(1), 'YLim',ylim);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function view = linearizeCmap(view)
%
% AUTHOR: Baseler/Poirson/Huk
% DATE:   3.4.99
% PURPOSE: Correct for the non linear input-output relationship
%	   between colormap entries and colors produced on the
%          the monitor.  This works best for 'hsv' colormaps.
% USAGE:
% HISTORY:  2.28.97 hab, poirson    Wrote mrLinearCmap
%           3.4.99 huk              Updated to mrLoadRet-2.0
% INPUT:   view (i.e. FLAT, INPLANE, or VOLUME)
% OUTPUT:  view with corrected colormap
%
% 2006/10 SD: updated so it works for every displayMode and colormap.

gammaLeft = 2;
gammaRight = 1.3;

cmap = eval(['view.ui.' view.ui.displayMode 'Mode.cmap']);
nGrays = eval(['view.ui.' view.ui.displayMode 'Mode.numGrays']);
nColors = eval(['view.ui.' view.ui.displayMode 'Mode.numColors']);

% different correction for different parts
iiLower = [nGrays+1:nGrays+round(nColors/2)];
iiUpper = [round(nGrays+nColors/2)+1:nGrays+nColors];

cmap(iiLower,:) = cmap(iiLower,:).^(1/gammaLeft);
cmap(iiUpper,:) = cmap(iiUpper,:).^(1/gammaRight);

colormap(cmap);

eval(['view.ui.' view.ui.displayMode 'Mode.cmap=cmap;']);

return


text:
function  peri = perimeter(grains,varargin)
% calculates the perimeter of a grain with or without inclusions
%
% Syntax
%
%   grains.perimeter
%   perimter(grains)
%   perimter(grains,'withInclusions')
%
% Input
%  grains - @grain2d
%
% Output
%  peri - perimeter (in measurement units)
%
% See also
% grain2d/equivalentPerimeter

if check_option(varargin,'withInclusion')

  bnd = grains.boundary;
  grainId = bnd.grainId;
  segLength = bnd.segLength;

  isGrain1 = grainId(:,1)>0;
  isGrain2 = grainId(:,2)>0;
  peri = accumarray(grainId(isGrain1,1),segLength(isGrain1),[max(grainId(:)),1]) ...
    + accumarray(grainId(isGrain2,2),segLength(isGrain2),[max(grainId(:)),1]);

  peri = peri(grains.id);

else

  poly = grains.poly;

  % remove inclusions
  incl = grains.inclusionId;
  for i = find(incl>0).'
    poly{i} = poly{i}(1:end-incl(i));
  end

  V = grains.V;

  peri =  cellfun(@(ind) sum(sqrt(sum(diff(V(ind,:)).^2,2))),poly);

end

author: mtex-toolbox

repo: mtex

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