text:
function linplus_test24 ( )

%*****************************************************************************80
%
%% TEST24 tests R8GB_ML.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 10;
  n = m;
  ml = 1;
  mu = 2;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST24\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  R8GB_ML computes A*x or A''*X\n' );
  fprintf ( 1, '    where A has been factored by R8GB_FA.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M              = %d\n', m );
  fprintf ( 1, '  Matrix columns N           = %d\n', n );
  fprintf ( 1, '  Lower bandwidth ML         = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU         = %d\n', mu );

  for job = 0 : 1
%
%  Set the matrix.
%
    [ a, seed ] = r8gb_random ( m, n, ml, mu, seed );
%
%  Set the desired solution.
%
    x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
    if ( job == 0 )
      b = r8gb_mxv ( m, n, ml, mu, a, x );
    else
      b = r8gb_vxm ( m, n, ml, mu, a, x );
    end
%
%  Factor the matrix.
%
    [ a_lu, pivot, info ] = r8gb_fa ( n, ml, mu, a );

    if ( info ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TEST24 - Fatal error!\n' );
      fprintf ( 1, '  R8GB_FA declares the matrix is singular!\n' );
      fprintf ( 1, '  The value of INFO is %d\n', info );
      return
    end
%
%  Now multiply factored matrix times solution to get right hand side again.
%
    b2 = r8gb_ml ( n, ml, mu, a_lu, pivot, x, job );

    if ( job == 0 )
      r8vec2_print_some ( n, b, b2, 10, '  A*x and PLU*x' );
    else
      r8vec2_print_some ( n, b, b2, 10, '  A''*x and (PLU)''*x' );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef Gumbel2C
%%GUMBEL2C A collection of functions for the bivariate Gumbel copulae.
%         Implemented functions include: PDF, CDF, tau.
%
%REFERENCES:
%[1] H. Joe and D. Kurowicka, Dependence Modeling: Vine Copula Handbook.
%    World Scientific, 2011.
%[2] C. Brechmann and U. Schepsmeier, "Dependence modeling with c- and
%    d-vine copulas: The r-package cdvine," Jan. 2011.
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    methods(Static)
        function jprob = PDF(u,theta,rot)
            %PDF Generates a joint probability for column vector u or
            %    a row vector of joint probabilities for matrix u
            %    with each column as the input vector.
            %
            %INPUTS:
            % u: A 2-by-n matrix where each column represents a
            %    2-dimensional random vector and n is the number of
            %    random vectors. Entries must be in [0,1].
            % theta: A scalar coupling parameter in the interval
            %        [1,inf).
            % rot: A scalar specifying one of four rotations of the
            %      standard bivariate Gumbel copula. Allowed values are 0,
            %      90, 180, and 270. Default is 0.
            %
            %OUTPUTS: jprob: A 1-by-n vector of joint probabilities.
            %
            %EXAMPLE 1: Make a contour plot of the density.
            %x = linspace(0,1,1e2);
            %[X,Y] = meshgrid(x);
            %Z = zeros(length(X),length(Y));
            %for idx = 1:length(X)
            %    for iidx = 1:length(Y)
            %        Z(idx,iidx) = Gumbel2C.PDF([X(idx,iidx);Y(idx,iidx)],3);
            %    end
            %end
            %contourf(X,Y,Z)
            %
            %EXAMPLE 2: Generate an exact joint density plot.
            % x = linspace(0,1,1e2);
            % [X,Y] = meshgrid(x);
            % Z = zeros(length(X),length(Y));
            % for idx = 1:length(X)
            % for iidx = 1:length(Y)
            % Z(idx,iidx) = Gumbel2C.PDF([X(idx,iidx);Y(idx,iidx)],3);
            % X(idx,iidx) = GumbelD.invCDF(X(idx,iidx),1,2);
            % Y(idx,iidx) = GaussianD.invCDF(Y(idx,iidx));
            % Z(idx,iidx) = Z(idx,iidx).*GumbelD.PDF(X(idx,iidx),1,2).*GaussianD.PDF(Y(idx,iidx));
            % end
            % end
            % surf(X,Y,Z)
            %
            %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
            %
            if ~exist('u','var') || isempty(u)
                jprob = [];
                return
            end
            if ~exist('theta','var') || isempty(theta)
                theta = 1;
            elseif theta < 1
                error('Theta parameter cannot be less than 1.')
            end
            if ~exist('rot','var') || isempty(rot)
                rot = 0;
            end

            if rot==0
                %No modifications necessary.
            elseif rot==90
                u(1,:) = 1-u(1,:);
            elseif rot==180
                u(1,:) = 1-u(1,:);
                u(2,:) = 1-u(2,:);
            elseif rot==270
                u(2,:) = 1-u(2,:);
            else
                error('Invalid rotation. Allowed values are 0, 90, 180, and 270.')
            end

            a = -log(u(1,:));
            b = -log(u(2,:));

            jprob = Gumbel2C.CDF(u,theta).*(a.*b).^(theta-1)...
                    .*((a.^theta+b.^theta).^(1/theta)+theta-1)...
                    ./(u(1,:).*u(2,:).*(a.^theta+b.^theta).^(2-(1/theta)));

            %NaN results from 0*inf in the calculation.
            %These should be 0s.
            nanIdxs = isnan(jprob);
            jprob(nanIdxs) = 0;

        end

        function jdist = CDF(u,theta,rot)
            %CDF Generates a joint distribution for column vector u or
            %    a row vector of joint probabilities for matrix u
            %    with each column as the input vector.
            %
            %INPUTS:
            % u: A 2-by-n matrix where each column represents a
            %    2-dimensional random vector and n is the number of
            %    random vectors. Entries must be in [0,1].
            % theta: A scalar coupling parameter in the interval
            %        [1,inf).
            % rot: A scalar specifying one of four rotations of the
            %      standard bivariate Gumbel copula. Allowed values are 0,
            %      90, 180, and 270. Default is 0.
            %
            %OUTPUTS: jprob: A 1-by-n vector of joint distribution values.
            %
            %EXAMPLE 1: Make a contour plot of the distribution.
            %x = linspace(0,1,1e2);
            %[X,Y] = meshgrid(x);
            %Z = zeros(length(X),length(Y));
            %for idx = 1:length(X)
            %    for iidx = 1:length(Y)
            %        Z(idx,iidx) = Gumbel2C.CDF([X(idx,iidx);Y(idx,iidx)],3);
            %    end
            %end
            %contourf(X,Y,Z)
            %
            %EXAMPLE 2: Generate an exact joint distribution plot.
            % x = linspace(0,1,1e2);
            % [X,Y] = meshgrid(x);
            % Z = zeros(length(X),length(Y));
            % for idx = 1:length(X)
            % for iidx = 1:length(Y)
            % Z(idx,iidx) = Gumbel2C.CDF([X(idx,iidx);Y(idx,iidx)],3);
            % X(idx,iidx) = GumbelD.invCDF(X(idx,iidx),1,2);
            % Y(idx,iidx) = GaussianD.invCDF(Y(idx,iidx));
            % Z(idx,iidx) = Z(idx,iidx);
            % end
            % end
            % surf(X,Y,Z)
            %
            %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
            %
            if ~exist('u','var') || isempty(u)
                jdist = [];
                return
            end
            if ~exist('theta','var') || isempty(theta)
                theta = 1;
            elseif theta < 1
                error('Theta parameter cannot be less than 1.')
            end
            if ~exist('rot','var') || isempty(rot)
                rot = 0;
            end

            if rot==0
                %No modifications necessary.
                a = -log(u(1,:));
                b = -log(u(2,:));

                jdist = exp(-(a.^theta+b.^theta).^(1/theta));
            elseif rot==90
                a = -log(1-u(1,:));
                b = -log(u(2,:));

                jdist = u(2,:)-exp(-(a.^theta+b.^theta).^(1/theta));
            elseif rot==180
                a = -log(1-u(1,:));
                b = -log(1-u(2,:));

                jdist = u(1,:)+u(2,:)-1+exp(-(a.^theta+b.^theta).^(1/theta));
            elseif rot==270
                a = -log(u(1,:));
                b = -log(1-u(2,:));

                jdist = u(1,:)-exp(-(a.^theta+b.^theta).^(1/theta));
            else
                error('Invalid rotation. Allowed values are 0, 90, 180, and 270.')
            end

        end

        function t = tau(theta)
            %%TAU Compute Kendall's tau for the copula.
            %
            %INPUTS:
            % theta: A scalar coupling parameter in the interval
            %        [1,inf).
            %
            %OUTPUTS:
            % t: The scalar value of Kendall's tau.
            %
            %October 2020 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
            %
            if ~exist('theta','var') || isempty(theta)
                theta = 1;
            elseif theta < 1
                error('Theta parameter cannot be less than 1.')
            end

            t = 1-1/theta;
        end
    end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function irepistruct = irepitiming(irepistruct)
% irepistruct = irepitiming(irepistruct)
%
% See irepistructure.m for more info on the struct parameters
%
% $Id: irepitiming.m,v 1.2 2016/04/11 03:19:08 greve Exp $

InvFlip = 180;

s = irepistruct; % copy into s for easy handling
s.tEvent = [];
s.EventSliceNo = []; % spatial/anatomical slice no
s.AcqSliceNo = []; % after inversion
s.FlipEvent = [];
s.IsReadOut = [];

t = 0;

% Dummies.
SliceOrderD = [1:s.nslices]; % Dummies do not permute
for nthtp = 1:s.ndummies
  if(s.Slice1PreInv==0)
    % Inversion
    t = t - s.PreInv;
    t = t + s.InvDur; % duration of inversion pulse
    tInv = t;
    s.tEvent = [s.tEvent; tInv];
    s.FlipEvent = [s.FlipEvent; InvFlip];
    s.IsReadOut = [s.IsReadOut; 0];
    s.TI = [s.TI; 0];
    s.AcqSliceNo = [s.AcqSliceNo; 0];
    s.EventSliceNo = [s.EventSliceNo; -1]; % Applies to all slices
    t = t + s.TI1;
  end
  % Readout (but no actual acq because these are dummies)
  nthslice = 0;
  for sliceno = SliceOrderD
    nthslice = nthslice + 1;
    s.tEvent = [s.tEvent; t];
    s.EventSliceNo = [s.EventSliceNo; sliceno];
    s.FlipEvent = [s.FlipEvent; s.ROFlip];
    s.IsReadOut = [s.IsReadOut; 0]; % 0 = no acq
    s.TI = [s.TI; 0];
    s.AcqSliceNo = [s.AcqSliceNo; 0];
    t = t + s.TBS;
    if(nthslice == 1 & s.Slice1PreInv)
      % Inversion
      t = t - s.PreInv;
      t = t + s.InvDur; % duration of inversion pulse
      tInv = t;
      s.tEvent = [s.tEvent; tInv];
      s.FlipEvent = [s.FlipEvent; InvFlip];
      s.IsReadOut = [s.IsReadOut; 0];
      s.TI = [s.TI; 0];
      s.AcqSliceNo = [s.AcqSliceNo; 0];
      s.EventSliceNo = [s.EventSliceNo; -1]; % Applies to all slices
      t = t + s.TI1;
    end
  end
end

% Non-dummies
for nthtp = 1:s.ntp
  if(s.Slice1PreInv==0)
    % Inversion before all slices
    t = t - s.PreInv;
    t = t + s.InvDur; % duration of inversion pulse
    tInv = t;
    s.tEvent = [s.tEvent; tInv];
    s.EventSliceNo = [s.EventSliceNo; -1]; % Applies to all slices
    s.FlipEvent = [s.FlipEvent; InvFlip];
    s.IsReadOut = [s.IsReadOut; 0];
    s.TI = [s.TI; 0];
    s.AcqSliceNo = [s.AcqSliceNo; 0];
    t = t + s.TI1;
  end
  % Read out
  SliceOrderR = circshift([1:s.nslices]',-(nthtp-1)*s.skip)';
  nthslice = 0;
  for sliceno = SliceOrderR
    nthslice = nthslice + 1;
    s.tEvent = [s.tEvent; t];
    s.EventSliceNo = [s.EventSliceNo; sliceno];
    s.FlipEvent = [s.FlipEvent; s.ROFlip];
    s.IsReadOut = [s.IsReadOut; 1];
    s.TI = [s.TI; t-tInv];
    s.AcqSliceNo = [s.AcqSliceNo; nthslice];
    t = t + s.TBS;
    if(nthslice == 1 & s.Slice1PreInv)
      % Inversion, when 1st slice occurs before inversion
      t = t - s.PreInv;
      t = t + s.InvDur;
      tInv = t;
      s.tEvent = [s.tEvent; tInv];
      s.EventSliceNo = [s.EventSliceNo; -1]; % Applies to all slices
      s.FlipEvent = [s.FlipEvent; InvFlip];
      s.IsReadOut = [s.IsReadOut; 0];
      s.TI = [s.TI; 0];
      s.AcqSliceNo = [s.AcqSliceNo; 0];
      t = t + s.TI1;
    end
  end % slice
end

irepistruct = s;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ f_rec , X ] = abel_inversion(h,R,upf,plot_results)
% This function calculates a Fourier-based Abel inversion based on the
% method described in [1]. Assuming cylindrical symmetry, the distribution
% function f_rec can be reconstructed from the measured profile h .
%
% The main idea is a Fourier-series-like expansion of the unknown
% distribution function f(r) into
% f(r) = sum_{n=lof}^{upf} (A_n * f_n(r))                               (1)
% where the lower frequency is set to 1 and the upper frequency upf is
% important for noise-filtering.
%
%
% INPUT:  H - dataset to be analyzed (1xN - matrix).
%           For optimum results, H(1) should be the center and H(N) the
%           edge of the investigated object and H(N) should be approx zero.
%           If no input H is given, a sampla dataset will be created.
%         R - radius of given system (i.e. distance between H(1) and H(N))
%         UPF  - upper frequency limit. Defines the number of cosinus
%           expansions (and critically determines computation time).
%           Choosing a very low value (e.g. 4) results in a low-pass
%           filtering effect, reducing noise (but also potential features)
%         PLOT_RESULTS - set to 1 to plot the results
%
% OUTPUT: F_REC - reconstructed density profile
%         X - x-vector containing spatial coordinates for F_REC
%
%
%     [1] G. Pretzler, Z. Naturforsch. 46a, 639 (1991)
%
% See also: COMPUTE_EXPANSION, GENERATE_TEST_DATA, SOLVE_LSQ
%
%                                         written by C. Killer, Sept. 2013

%% format data / generate sample data

% create sample data based on a polynomial distribution function
% if no data input is given
if ~exist('h', 'var') || isempty(h)
    [X,h,R]=generate_test_data;
    plot_results=1;
else
    X=linspace(0,R-0.01,length(h))';
end

% default value for number of expansion elements
if ~exist('upf', 'var'); upf=10; end;

% avoiding problems if this flag is not given in input
if ~exist('plot_results', 'var'); plot_results=0; end;

%% calculate series expansions fn and corresponding integrals hn

[fn,hn] = compute_expansion( X,upf,R );

%% Solving optimization problem

% guess some initial values for optimisation
x0=1*ones(upf+1,1);

% solve for amplitudes A
A = solve_lsq(h,hn,x0);

% create vector for resulting reconstructed distribution
f_rec=zeros(length(h),1);

% special case for n=0 (where f_0(r) = 1)
f_rec = f_rec + A(1)*1;

% iterate eq. (1) for n=1:upf
for c=2:length(x0)
    f_rec = f_rec + A(c).*fn(:,c);
end

if plot_results==1
    figure; % normalized profiles for better comparison
    set(gca,'linewidth',1.5,'fontsize',16)
    plot(X,f_rec./max(f_rec),'k','Linewidth',1.5);
    hold on; plot(X,h./max(h),'b','Linewidth',1.5);
    grid on; box on;
    title(sprintf('number of cos-expansions: %i',upf))
    legend('reconstructed distribution','measured profile','Location','SouthWest')
end

author: Sable

repo: mcbench-benchmarks


text:
function wfun = estimateWBTVPriorWeights

    wfun.function = @WBTVPriorWeightsWeightingFunction;
    wfun.parameters = {0.5, ... % Sparsity parameter
                       2};      % Tuning constant

function [weights, scaleParameter] = WBTVPriorWeightsWeightingFunction(z, weights, sparsityParam, tuningConstant)

    if nargin < 3
        % Use default value p = 0.5 for the sparsity parameter.
        sparsityParam = 0.5;
    end
    if nargin < 4
        tuningConstant = 2;
    end

    % Adaptive estimation of the scale parameter.
    scaleParameter = getAdaptiveScaleParameter(z, weights);

    % Estimation of the weights based on pre-selected scale parameter.
    weights = (sparsityParam * (tuningConstant*scaleParameter)^(1-sparsityParam)) ./ (abs(z).^(1-sparsityParam));
    weights(abs(z) <= (tuningConstant*scaleParameter)) = 1;

function scaleParameter = getAdaptiveScaleParameter(z, weights)

    scaleParameter = weightedMedian( abs(z - weightedMedian(z, weights)), weights );

author: thomas-koehler

repo: SupER


text:
function [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval )

%*****************************************************************************80
%
%% I4VEC_INDEX_INSERT inserts a value in an indexed sorted I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the size of the current list.
%
%    Input, integer X(N), the list.
%
%    Input, integer INDX(N), the sort index of the list.
%
%    Input, integer XVAL, the value to be sought.
%
%    Output, integer N, the size of the current list.
%
%    Output, integer X(N), the list.
%
%    Output, integer INDX(N), the sort index of the list.
%
  if ( n <= 0 )
    n = 1;
    x = [ xval ];
    indx = [ 1 ];
    return
  end

  [ less, equal, more ] = i4vec_index_search ( n, x, indx, xval );

  x(n+1) = xval;
  indx(n+1:-1:more+1) = indx(n:-1:more);
  indx(more) = n + 1;
  n = n + 1;

  return
end


text:
function [ vX ] = SolveL1NormSetMinimization( mY )
% ----------------------------------------------------------------------------------------------- %
% [ vX, mX ] = SolveLsBoxConstraints( mA, vB, vC, vD, vX, numIterations, stopTol )
%   Solves \arg \min_{x} 0.5 || A x - b ||, s.t. c <= x <= d using
%   Gradient Descent Method.
% Input:
%   - mA            -   Model Matrix.
%                       Input model matrix.
%                       Structure: Matrix (m x n).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% Output:
%   - vX            -   Solution Vector.
%                       The solution to the optimization problem..
%                       Structure: Vector (m x 1).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  A
% Remarks:
%   1.  B
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     02/03/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

vecDim      = size(mY, 1);
numSamples  = size(mY, 2);

% vX = [vX; vT]
numElements = vecDim * numSamples;

% Summing over all {t}_{ij}
vF = [zeros(vecDim, 1); ones(numElements, 1)];

% Building the matrix where we run on {x}_{i}, {t}_{i, j}, {y}_{i, j} in
% column wise manner. Where i is the element index and j is the sample
% index.
mA = [repmat(speye(vecDim), numSamples, 1), -speye(numSamples * vecDim); repmat(-speye(vecDim), numSamples, 1), -speye(numSamples * vecDim)];
vB = [mY(:); -mY(:)];

sSolverOptions = optimoptions('linprog', 'Display', 'off');
vX = linprog(vF, mA, vB, [], [], [], [], sSolverOptions);
vX = vX(1:vecDim);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
% -----------------------------------------------------------------------------------------------------
function pyramid = make_scale_pyramid(im, targetPosition, in_side_scaled, out_side, avgChans, stats, p)
%MAKE_SCALE_PYRAMID
%   computes a pyramid of re-scaled copies of the target (centered on TARGETPOSITION)
%   and resizes them to OUT_SIDE. If crops exceed image boundaries they are padded with AVGCHANS.
%
%   Luca Bertinetto, Jack Valmadre, Joao F. Henriques, 2016
% -----------------------------------------------------------------------------------------------------
    in_side_scaled = round(in_side_scaled);
    pyramid = gpuArray(zeros(out_side, out_side, 3, p.numScale, 'single'));
    max_target_side = in_side_scaled(end);
    min_target_side = in_side_scaled(1);
    beta = out_side / min_target_side;
    % size_in_search_area = beta * size_in_image
    % e.g. out_side = beta * min_target_side
    search_side = round(beta * max_target_side);
    [search_region, ~] = get_subwindow_tracking(im, targetPosition, [search_side search_side], [max_target_side max_target_side], avgChans);
    if p.subMean
        search_region = bsxfun(@minus, search_region, reshape(stats.x.rgbMean, [1 1 3]));
    end
    assert(round(beta * min_target_side)==out_side);

    for s = 1:p.numScale
        target_side = round(beta * in_side_scaled(s));
        pyramid(:,:,:,s) = get_subwindow_tracking(search_region, (1+search_side*[1 1])/2, [out_side out_side], target_side*[1 1], avgChans);
    end
end

author: bertinetto

repo: siamese-fc


text:
function cvx_expert_check( fname, varargin )

global cvx___
if cvx___.expert, return; end
if ~isempty( varargin ) && ~any(cellfun('isclass',varargin,'cvx')), return; end

url = [ 'file:///', cvx___.where, cvx___.fs, 'doc', cvx___.fs, 'advanced.html#the-successive-approximation-method' ];
fprintf( 1, [ 'CVX Warning:\n', ...
'   Models involving "%s" or other functions in the log, exp, and entropy\n', ...
'   family are solved using an experimental successive approximation method.\n', ...
'   This method is slower and less reliable than the method CVX employs for\n', ...
'   other models. Please see the section of the user''s guide entitled\n', ...
'       The successive approximation method\n', ...
'   for more details about the approach, and for instructions on how to\n', ...
'   suppress this warning message in the future.\n' ], fname, url );

cvx___.expert = true;

% Note that we do *not* call cvx_save_prefs here. We only save the
% preferences if an explicit setting of "cvx_expert true" is made.

% 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 [hPa] = Pa2hPa(Pa)
% Convert pressure from pascals to hectopascals.
% Chad Greene 2012
hPa = Pa*0.010;


text:
function circle_segment_test11 ( )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_TEST11 demonstrates CIRCLE_SEGMENT_ROTATION_FROM_CHORD.
%
%  Discussion:
%
%    We make a table of all pairs of angles that are multiples of pi/12.
%
%    For each pair, we compute the rotation, that is, the angle of the
%    central radius of the circle segment.  We print out the result in
%    terms of multiples of pi/12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_SEGMENT_TEST11:\n' );
  fprintf ( 1, '  CIRCLE_SEGMENT_ROTATION_FROM_CHORD is given the endpoints\n' );
  fprintf ( 1, '  of a chord, and is asked to determine the angle of the\n' );
  fprintf ( 1, '  central radius vector.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We make a table of all pairs of angles that are multiples\n' );
  fprintf ( 1, '  of pi/12, determine the corresponding chord endpoints, and\n' );
  fprintf ( 1, '  compute the rotation angle, also printed as a multiple of pi/12.\n' );

  r = 2.0;
  c = [ 3.0; 5.0 ];
  fprintf ( 1, '\n' );
  fprintf ( 1, '            0.0   1.0   2.0   3.0   4.0   5.0   6.0   7.0   8.0   9.0  10.0  11.0  12.0\n' );
  fprintf ( 1, '\n' );
  for i = 0 : 12
    rho1 = i * pi / 6;
    p1 = c + r * [ cos(rho1); sin(rho1) ];
    fprintf ( 1, '  %2d.0:  ', i );
    for j = 0 : 12
      rho2 = j * pi / 6;
      p2 = c + r * [ cos(rho2); sin(rho2) ];
      alpha = circle_segment_rotation_from_chord ( r, c, p1, p2 );
      fprintf ( 1, '  %4.1f', round ( 10 * 6 * alpha / pi ) / 10 );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function dat=readinr(fname)
%
% vol=readinr(fname)
%
% load a volume from an INR file
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2009/05/03
%
% input:
%      fname: input file name
%
% output:
%      dat: output, data read from the inr file
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fid=fopen(fname,'rb');
s=fread(fid,256,'uchar');

s=char(s)';

if(regexp(s,'#INRIMAGE-4')~=1)
   error('INRIMAGE header was not found')
end

nx=regexp(s,'XDIM\s*=\s*([0-9]+)','tokens');
if(length(nx))
   nx=str2num(nx{1}{1});
else
   error('no XDIM found');
end

ny=regexp(s,'YDIM\s*=\s*([0-9]+)','tokens');
if(length(ny))
   ny=str2num(ny{1}{1});
else
   error('no YDIM found');
end

nz=regexp(s,'ZDIM\s*=\s*([0-9]+)','tokens');
if(length(nz))
   nz=str2num(nz{1}{1});
else
   error('no ZDIM found');
end

nv=regexp(s,'VDIM\s*=\s*([0-9]+)','tokens');
if(length(nv))
   nv=str2num(nv{1}{1});
else
   nv=1;
end

type=regexp(s,'TYPE=([a-z ]+)','tokens');
if(length(type))
   type=type{1}{1};
else
   error('no TYPE found');
end

pixel=regexp(s,'PIXSIZE=([0-9]+)','tokens');
if(length(pixel))
   pixel=str2num(pixel{1}{1});
else
   error('no PIXSIZE found');
end

%header=sprintf(['#INRIMAGE-4#{\nXDIM=%d\nYDIM=%d\nZDIM=%d\nVDIM=1\nTYPE=%s\n' ...
%  'PIXSIZE=%d bits\nCPU=decm\nVX=1\nVY=1\nVZ=1\n'],size(vol),btype,bitlen);

if(strcmp(type,'unsigned fixed') & pixel==8)
   dtype='uint8';
elseif(strcmp(type,'unsigned fixed') & pixel==16)
   dtype='uint16';
elseif(strcmp(type,'float') & pixel==32)
   dtype='float';
elseif(strcmp(type,'float') & pixel==64)
   dtype='double';
else
   error('volume format not supported');
end

dat=fread(fid,nx*ny*nz*nv,dtype);
fclose(fid);

if(nv==1)
   dat=reshape(dat,[nx,ny,nz]);
else
   dat=reshape(dat,[nx,ny,nz,nv]);
end

repo: REST


text:
%% A Simple and Robust Vertical Handoff Algorithm for HeterogeneousWireless Mobile Networks
%   Function Parameters:
%   C=cost of Ith network      S=Security     P=Power
%   D = Network Conditions     F= Network Performance
%   V = Velocity of MN         T= Estimated Time MN will stay in the
%                                   network coverage
%   NS = Network Set of eligible candidates for VHO
%   N = Number of Networks Presently available for evaluation.
%   lambda = mean number of request arrivals per unit time
%   mu = mean number of calls serviced per unit time
%   wx = weight function of the 'x' parameter for evaluating EVHDF
%
% Take as input from user the various network parameters for each network
% under consideration.
%
%   Reference: A Simple and Robust Vertical Handoff Algorithm for
%   Heterogenous Wireless Mobile Networks
%   Paper Authors: Daojing He, Caixia Chi, Sammy Chan, Chun Chen, Jiajun Bu,
%   Mingjian Yin.
%   Wireless Pers Commun DOI 10.1007/s11277-010-9922-x
%
%   Code Author: Anshul Thakur
%   Contact: anshulthakurjourneyendless@gmail.com
%   Created: 4/4/2011
%   Copyright Owner: Anshul Thakur

%% Note that the code that has been commented out are of the various
% parameters that the author did not consider for his evaluation. User may
% uncomment the necessary part to include more parameters in their
% evaluation
function code_modified()
%This section of code takes in various parameters of the different networks
%between which decision of a handoff is to be made.
fprintf('\nThe number of networks under evaluation N: \n')
N=input('=');
b=zeros(1,N);
RSS=zeros(1,N);
%V=zeros(1,N);
T=zeros(1,N);
%P=zeros(1,N);
%C=zeros(1,N);
%fj=zeros(1,N);
%wc=zeros(1,N);
%ws=zeros(1,N);
%wf=zeros(1,N);
wp=zeros(1,N);
wd=zeros(1,N);
M=zeros(1,N);
lambda=zeros(1,N);
mu=zeros(1,N);
pj=zeros(1,N);
NL=0;
fprintf('\nEnter the threshold values of various parameters when prompted.\n');
    bi=input('Threshold Current Available Bandwidth: ');
    RSSi=input('Threshold Received Signal Strength:  ');
    %Vi=input('Threshold Velocity of Mobile Station:  ');
    Ti=input('Threshold Estimated Time MS will be in present network:  ');
    %Pi=input('Threshold Battery Power of MS  :');
    %Ci=input('Threshold Cost of Network to MS  :');
for i=1:N
    fprintf('Enter the values for the network number %d when prompted.\n', i);
    if(i==1)
        fprintf('Enter the values of the currently connected network\n');
    end
    b(1,i)=input('Current Available Bandwidth: ');
    RSS(1,i)=input('Received Signal Strength:  ');
    %V(1,i)=input('Velocity of Mobile Station:  ');
    T(1,i)=input('Estimated Time MS will be in present network:  ');
    %P(1,i)=input('Battery Power of MS  :');
    %C(1,i)=input('Cost of Network to MS  :');
    %s(1,i)=input('Security level in Network    :');
    pj(1,i)=input('Power Dissipation in Network  :');
    %fj(1,i)=input('Network Performance Parameter of Network    :');
    lambda(1,i)=input('Mean number of request arrivals per unit time    :');
    mu(1,i)=input('Mean number of calls serviced per unit time  :');
    fprintf('\n \t Enter the Network Dependent Weights to the following parameters\n');
    %wc(1,i)=input('Cost of Service:');
    %ws(1,i)=input('Security Parameter:');
    wp(1,i)=input('Power Consumption: ');
    wd(1,i)=input('Network Conditions:  ');
    %wf(1,i)=input('Network Performance     :');

    %% This part evaluates the first step decision function to find the eligible
    % set of networks out of the available networks with values greater
    % than threshold
    %M(1,i)=((b(1,i)-bi)>0)&&((RSS(1,i)-RSSi)>0)&&((V(1,i)-Vi)>0)&&((T(1,i)-Ti)>0)&&((P(1,i)-Pi)>0)&&((C(1,i)-Ci)>0);
    M(1,i)=((b(1,i)-bi)>0)&&((RSS(1,i)-RSSi)>0)&&((T(1,i)-Ti)>0);
    if M(1,i)==1
        NL=NL+1;
    end
end
if(NL==0)
    fprintf('\nNo eligible networks found. Remain Connected to existing network itself\n');
else
%% Decision making when Network list has some networks in it
    NS=zeros(1,NL);
    hi=zeros(1,NL);
    EQ=zeros(1,NL);
    Dj=zeros(1,NL);
   x=1;
   for i=1:N
        if M(1,i)~=0
            NS(1,x)=i;
            x=x+1;
        end
   end
   if ((NL==1)&&(NS(1,1)==1))
       fprintf('\nNo other network is having minimum service quality better than threshold. Stay in same network\n');
   elseif((NL==1)&&(NS(1,1)~=1))
       fprintf('\n Handoff to new network with Network ID %d',NS(1,1));
   else
       fprintf('\nIs the Mobile Node Reource rich or Resource Poor?\n');
       res=input('r for Rich/ p for Poor      ','s');
       if(strcmp(res,'p'))
           if(NS(1,1)==1)
               frpintf('\nNo handoff required for Resource Poor Mobile Node at present\n');
           else
               fprintf('\nPerform Handoff with any of the %d Networks in the pool',NL);
           end
       elseif(strcmp(res,'r'))
           for i=1:NL
               index=NS(1,i);
               hi(1,i)=dyn_cal_prob(lambda(1,index),mu(1,index),b(1,index));
               Dj(1,i)=(b(1,index))/(hi(1,i));
           end
           %% This section evaluates the Extended Vertical Handoff Decision
           %  Function for the network pool previously selected as eligible
           for i=1:NL
               %index=NS(i,i);
               %Ct(1,i)=C(1,index);
               %EQ(1,i)=((wc(1,index)*(1/Ct(1,i)))/min(Ct,[],2))+(ws(1,index)*s(1,index)/max(s,[],2))+(wp(1,index)*pj(1,index)/max(pj,[],2))+(wd(1,index)*Dj(1,i)/max(Dj,[],2))+(wf(1,index)*fj(1,index)/max(fj,[],2));
               %EQ(1,i)=Dj(1,i)/max(Dj,[],2);
               EQ(1,i)=(wp(1,index)*pj(1,index)/max(pj,[],2))+(wd(1,index)*Dj(1,i)/max(Dj,[],2));
           end
           [para,netwrk]=max(EQ,[],2);
           if netwrk==1
               fprintf('\nCurrently selected Network is the best network, no Handoff required\n');
           else
               fprintf('\nPerform handoff with Network %d for best performance\n',NS(netwrk));
           end
       end
   end
    figure;
    subplot(2,1,1);
    stem(NS(1,:),hi(1,:));
    ylabel('New Dynamic Call Blocking Probability');
    xlabel('Network Number');
    title('a.   NDCBP vs Network ');
    subplot(2,1,2);
    stem(NS(1,:),EQ(1,:));
    title('   b.  EVHDF vs Network');
    xlabel('Netowrk');
    ylabel('Extended Vertical Handoff Decision Function');
end
end

%% Function to compute the New Dynamic Call blocking Probability
function hi=dyn_cal_prob(lambda, mu, bi)
    b=(lambda/mu);
    a=(b.^(bi))/factorial(bi);
    sumi=0;
    for n=0:bi
        sumi=sumi+((b.^exp(n))/factorial(n));
    end
    hi=a/sumi;
end


text:
function [h,p,T] = BartlettTest(data,alpha)

%BartlettTest - Test if k groups of samples have equal variances (homogeneity of variances).
%
%  This test assumes that all populations follow a Gaussian distribution.
%
%  USAGE
%
%    [h,p,t] = BartlettTest(data,alpha)
%
%    data           Nx2 matrix of (observation,group) pairs
%    alpha          optional significance level (default = 0.05)
%
%    h              test result (1 = reject null hypothesis, 0 = accept)
%    p              p value
%    t              test statistics

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% 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.

% Significance level
if nargin < 2,
	alpha = 0.05;
end

% Number of groups
k = max(data(:,2));
for i = 1:k,
	group = data(:,2) == i;
	% Number of samples in this group
	n(i) = sum(group);
	% Variance for this group
	s2(i) = var(data(group,1));
end
% Total number of samples
N = sum(n);
% Pooled variance
S2 = sum((n-1).*s2)/(N-k);

% Test statistics (biased)
t = (N-k)*log(S2)-sum((n-1).*log(s2));
% Bias correction
C = 1+(1/(3*(k-1)))*(sum(1./(n-1))-1/(N-k));
% Test statistics (unbiased)
T = t/C;

% Chi square at alpha with (k-1) degrees of freedom
p = 1 - chi2cdf(T,k-1);

h = p < alpha;

disp(['Bartlett test: Variances: ' num2str(s2)]);

if h,
	message = '+++ Two or more variances are significantly different';
else
	message = '--- Variances are not significantly different';
end
message = ['Bartlett test: ' message ' (p='  num2str(p) ', T=' num2str(T) ', N=' int2str(N) ')'];
disp(message);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function b = cppsl ( ap, n, b )

%*****************************************************************************80
%
%% CPPSL solves a complex hermitian positive definite linear system.
%
%  Discussion:
%
%    The matrix is assumed to have been factored by CPPCO or CPPFA.
%
%    A division by zero will occur if the input factor contains
%    a zero on the diagonal.  Technically this indicates
%    singularity but it is usually caused by improper subroutine
%    arguments.  It will not occur if the subroutines are called
%    correctly and INFO == 0.
%
%    To compute inverse(A) * C where C is a matrix with P columns:
%
%      call cppco(ap,n,rcond,z,info)
%
%      if (rcond is too small .or. info /= 0) then
%        error
%      end if
%
%      do j = 1, p
%        call cppsl(ap,n,c(1,j))
%      end do
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex AP(N*(N+1)/2), the output from CPPCO or CPPFA.
%
%    Input, integer N, the order of the matrix.
%
%    Input, complex B(N), the right hand side.
%
%    Output, complex B(N), the solution.
%
  kk = 0;
  for k = 1 : n
    t = conj ( ap(kk+1:kk+k-1) ) * transpose ( b(1:k-1) );
    kk = kk + k;
    b(k) = ( b(k) - t ) / ap(kk);
  end

  for k = n : -1 : 1
    b(k) = b(k) / ap(kk);
    kk = kk - k;
    t = -b(k);
    b(1:k-1) = b(1:k-1) + t * ap(kk+1:kk+k-1);
  end

  return
end


text:
clear all
videofile = 'viptraffic.avi';
info = mmfileinfo(videofile);
cols=info.Video.Width;
rows=info.Video.Height;
hReader = vision.VideoFileReader(videofile,...
    'ImageColorSpace', 'RGB',...
    'VideoOutputDataType', 'single');
hFlow = vision.OpticalFlow( ...
    'OutputValue', 'Horizontal and vertical components in complex form', ...
    'ReferenceFrameDelay', 3,...
    'Method','Horn-Schunck');
hMean1 = vision.Mean;
hMean2 = vision.Mean('RunningMean', true);
hFilter = vision.MedianFilter;
hClose = vision.MorphologicalClose('Neighborhood', strel('line',5,45));
hBlob = vision.BlobAnalysis(...
    'CentroidOutputPort', false,...
    'AreaOutputPort', true, ...
    'BoundingBoxOutputPort', true,...
    'OutputDataType', 'double', ...
    'MinimumBlobArea', 250,...
    'MaximumBlobArea', 3600,...
    'MaximumCount', 80);
hErode = vision.MorphologicalErode('Neighborhood', strel('square',2));
hShape1 = vision.ShapeInserter(...
    'BorderColor', 'Custom', ...
    'CustomBorderColor', [0 1 0]);
hShape2 = vision.ShapeInserter(...
    'Shape','Lines', ...
    'BorderColor', 'Custom', ...
    'CustomBorderColor', [255 255 0]);
hText = vision.TextInserter(...
    'Text', '%4d',...
    'Location',  [1 1], ...
    'Color', [1 1 1],...
    'FontSize', 12);
sz = get(0,'ScreenSize');
pos = [(sz(3)-4*(cols+75))/2, (sz(4)-rows)/2 cols+60 rows+80];
hVideo1 = vision.VideoPlayer('Name','Original Video','Position',pos);
pos(1) = pos(1)+cols+75;
hVideo2 = vision.VideoPlayer('Name','Motion Vector','Position',pos);
pos(1) = pos(1)+cols+75;
hVideo3 = vision.VideoPlayer('Name','Thresholded Video','Position',pos);
pos(1) = pos(1)+cols+75;
hVideo4 = vision.VideoPlayer('Name','Results Video','Position',pos);
[xpos,ypos]=meshgrid(1:5:cols,1:5:rows);
xpos=xpos(:);
ypos=ypos(:);
locs=sub2ind([rows,cols],ypos,xpos);
while ~isDone(hReader)
    pause(0.3);
    frame  = step(hReader);
    gray = rgb2gray(frame);
    flow = step(hFlow, gray);
    lines = [xpos, ypos, xpos+20*real(flow(locs)), ypos+20*imag(flow(locs))];
    vector = step(hShape2, frame, lines);
    magnitude = flow .* conj(flow);
    threshold = 0.5 * step(hMean2, step(hMean1, magnitude));
    carobj = step(hFilter, magnitude >= threshold);
    carobj = step(hClose, step(hErode, carobj));
    [area, bbox] = step(hBlob, carobj);
    grow=22;
    idx = bbox(:,1) > grow;
    ratio = zeros(length(idx), 1);
    ratio(idx) = single(area(idx,1))./single(bbox(idx,3).*bbox(idx,4));
    flag = ratio > 0.4;
    count = int32(sum(flag));
    bbox(~flag, :) = int32(-1);
    result = step(hShape1, frame, bbox);
    result(grow:grow+1,:,:) = 1;
    result(1:15,1:30,:) = 0;
    result = step(hText, result, count);
    step(hVideo1, frame);
    step(hVideo2, vector);
    step(hVideo3, carobj);
    step(hVideo4, result);
end
release(hReader);

author: decouples


text:
function [d,BMA,PEBs] = spm_dcm_bdc(GCMs,field,M,ynames,models,noplot)
% Compare datasets using DCM and PEB (Bayesian data comparison)
%
% Performs the following procedure:
%
% 1. Identifies the optimum reduced model (PEB) collapsed over datasets.
% 2. For each dataset, creates reduced GCMs based on the optimum model
%    above and estimates a per-dataset PEB model.
% 3. Computes a series of measures on each dataset's PEB model
%
% Inputs:
%
% GCMs   - {1 x ny} cell array of filenames, each of which is a GCM file.
% field  - {1 x nf} cell array of field names, e.g. {'A','B'}
% M      - (Optional) struct for configuring PEB. See spm_dcm_peb.m .
% M.gamma- (Optional) prior variance for switched off parameters when
%          automatically forming a model space for KL-over-models.
% M.bmr  - (Optional) if true, searches for an optimal model across all
%          datasets before performing BDC. If false, all parameters from
%          the DCM are included. [default: true];
% ynames - (Optional) {1 x ny} cell array of names for each dataset
% models - (Optional) model space to use for computing the information gain
%          over models. Accepts a nested cell array of parameter names, or
%          a binary matrix (models x parameters) with which parameters
%          to switch on or off in each model.
% noplot - (Optional) if true, does not show plots.
%
% ny = number of datasets, nf = number of DCM fields, nm = number of models
%
% Returns:
%
% d.precisions - [np x ny] Precision of each DCM parameter in each dataset
% d.dcm_negent - [1 x ny]  Negative entropy (certainty) of DCM parameters
% d.rfx_negent - [1 x ny]  Negative entropy (certainty) of the estimated
%                          between-subject variability
% d.complexity - [1 x ny]  Number of effective parameters in the model
% d.model_F    - [nm x ny] Free energy of each candidate model
% d.model_P    - [nm x ny] Posterior probability of each candidate model
% d.model_KL   - [1 x ny]  Ability to disriminate between similar models
%
% BMA        - Bayesian model average across all datasets
% PEBs       - [1 x ny] PEB for each dataset
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

% Peter Zeidman & Samira Kazan
% $Id: spm_dcm_bdc.m 7495 2018-11-22 15:52:07Z peter $

if nargin < 5
    models = [];
end

if nargin < 6
    noplot = false;
end

% Load and concatenate models across datasets
% -------------------------------------------------------------------------
GCM = {};
for i = 1:length(GCMs)
    if ischar(GCMs{i})
        gcm = load(GCMs{i});
        gcm = gcm.GCM;
    else
        gcm = GCMs{i};
    end

    GCM = [GCM; gcm(:,1)];
end

ns = size(gcm,1);   % Number of subjects
ny = length(GCMs);  % Number of datasets
nm = ns*ny;         % Number of DCMs

% Set defaults
% -------------------------------------------------------------------------
if nargin < 2 || isempty(field), field = {'B'}; end
if nargin < 3, M = []; end
if nargin < 4 || isempty(ynames)
    ynames = cell(1,ny);
    for i = 1:ny
        ynames{i} = sprintf('Dataset %d',i);
    end
end

if isempty(M)
    M       = struct();
    M.alpha = 1;
    M.beta  = 16;
    M.hE    = 0;
    M.hC    = 1/16;
    M.Q     = 'single';
    M.X     = ones(nm,1);
end

try
    bmr = M.bmr;
catch
    bmr = true;
end

try
    gamma = M.gamma;
catch
    gamma = 0; % Point null hypothesis
end

% Identify optimum model structure across all datasets
% -------------------------------------------------------------------------
if bmr
    % Perform PEB analysis
    PEB = spm_dcm_peb(GCM,M,field);

    % Prune full model
    BMA = spm_dcm_peb_bmc(PEB);

    % Identify disabled PEB parameters
    is_disabled = find(BMA.Pp < 1/1000);

    % Get indices of disabled parameters
    disabled_pidx = BMA.Pind(is_disabled);
else
    disabled_pidx = [];
    BMA = [];
end

% Estimate PEBs for each dataset using the optimal DCM structure from above
% -------------------------------------------------------------------------

% Update design matrix
M.X = ones(ns,1);

PEBs = {};
for i = 1:ny
    % Load dataset
    GCM = spm_dcm_load(GCMs{i});

    % Disable any parameters which were pruned from the PEB.
    rE = GCM{1}.M.pE;
    rC = GCM{1}.M.pC;
    for j = 1:length(disabled_pidx)
        rC(disabled_pidx(j),disabled_pidx(j)) = 0;
    end

    % Apply priors at DCM level
    GCM = spm_dcm_reduce(GCM, rE, rC);

    % Build PEB
    M.X = ones(size(GCM,1),1);
    PEBs{i} = spm_dcm_peb(GCM,M,field);

end

% Compute measures
% -------------------------------------------------------------------------
d = struct();
d.precisions = [];
d.dcm_negent = [];
d.rfx_negent = [];
d.complexity = [];

for i = 1:ny
    % Certainty about RFX
    d.rfx_negent(i) = -0.5 * spm_logdet(2 * pi * exp(1) * PEBs{i}.Ch);

    % Certainty about parameters
    d.dcm_negent(i) = -0.5 * spm_logdet(2 * pi * exp(1) * PEBs{i}.Cp);

    % Information gain (KL-divergence) over parameters
    qE    = PEBs{i}.Ep;    % Posterior expectations
    pE    = PEBs{i}.M.pE;  % Prior expectations
    qC    = PEBs{i}.Cp;    % Posterior covariance
    pC    = PEBs{i}.M.pC;  % Prior covariance

    d.complexity(i) = spm_kl_normal(qE,qC,pE,pC);
end

% Assess model discrimination
% -------------------------------------------------------------------------

% If no models are provided, enter an automatic mode where difficult to
% distinguish parameters are identified
is_automatic = isempty(models);

% Form a model space K with one row per model and one column per parameter
% where 1 = on and 0 = off.
if isempty(models)
    % Default: each parameter turned off individually
    np = length(PEBs{1}.Pnames);
    K  = ones(np, np);
    for i = 1:np
        K(i,i) = 0;
    end
elseif iscell(models) && ~isempty(models) && iscell(models{1})
    % Nested cell array of included fields
    Pnames = char(PEBs{1}.Pnames);

    K  = [];

    for m = 1:length(models)
        k = any(ismember(Pnames,models{m}),2);
        K(end+1,:) = k';
    end

    % Remove any models where parameters were provided but none was found
    % in the PEB.
    to_remove = ~cellfun(@isempty,models) & ~any(K,2);
    K(to_remove,:) = [];
elseif ismatrix(models)
    % Pre-defined matrix of parameters
    K = models;
else
    error('Unknown model specification');
end

% Calculate evidence for reduced models with each connection switched off
% Fs(m,i) is the free energy of the PEB model with parameter m disabled
% in dataset i.
Fs = []; % models x datasets
for i = 1:ny
    qE = PEBs{i}.Ep;
    qC = PEBs{i}.Cp;
    pE = PEBs{i}.M.pE;
    pC = PEBs{i}.M.pC;

    for m = 1:size(K,1)
        rE = pE;
        rC = pC;
        off = find(K(m,:) == 0);
        rC(off,off) = gamma;

        Fs(m,i) = spm_log_evidence_reduce(qE,qC,pE,pC,rE,rC);
    end
end

if is_automatic
    % Retain nr difficult to discriminate models
    i  = find(mean(Fs,2) > -3);
else
    % Retain them all
    i = 1:size(Fs,1);
end
nr = numel(i);

% Evidence -> probability
p  = spm_softmax(Fs(i,:));

% KL divergence (information gain)
d.model_KL = sum(p.*log(p)) + log(nr);

d.model_F = Fs;
d.model_P = p;

% Get PEB parameter precisions for each dataset
% -------------------------------------------------------------------------
for i = 1:ny
    d.precisions(:,i) = diag(spm_inv(PEBs{i}.Cp));
end

if noplot
    return;
end

% Plot PEB parameters for each dataset
% -------------------------------------------------------------------------
spm_figure('GetWin','Datasets');
rows = max(ceil(ny / 2),3);
cols = 2;

for i = 1:ny
    subplot(rows,cols,i);
    spm_plot_ci(PEBs{i}.Ep,PEBs{i}.Cp);
    title(ynames{i},'FontSize',16);
    xlabel('PEB Parameter','FontSize',12);
    ylabel('Estimate','FontSize',12);
    ylim([-2 2]);
end

% Plot measures
% -------------------------------------------------------------------------
spm_figure('GetWin','Data Comparison');
spm_clf;
rows = 3; cols = 2;

subplot(rows,cols,1:2);
bar(d.precisions);
title('Parameter certainty','FontSize',16);
xlabel('Parameter','FontSize',12); ylabel('Precision','FontSize',12);
legend(ynames,'Location','South','Orientation','horizontal','FontSize',12);

subplot(rows,cols,3);
y = d.dcm_negent - min(d.dcm_negent);
bar(y);
hold on; p = scatter(1:length(y),y); set(p,'visible','off'); hold off
title('Parameter certainty','FontSize',16);
xlabel('Dataset','FontSize',12); ylabel('Relative neg entropy (nats)','FontSize',12);

subplot(rows,cols,4);
y = d.rfx_negent - min(d.rfx_negent);
bar(y);
hold on; p = scatter(1:length(y),y); set(p,'visible','off'); hold off
title('RFX certainty','FontSize',16);
xlabel('Dataset','FontSize',12); ylabel('Relative neg entropy (nats)','FontSize',12);

subplot(rows,cols,5);
y = d.complexity-min(d.complexity);
bar(y);
hold on; p = scatter(1:length(y),y); set(p,'visible','off'); hold off
title('Information gain (parameters)','FontSize',14);
xlabel('Dataset','FontSize',12); ylabel('Relative neg entropy (nats)','FontSize',12);

subplot(rows,cols,6);
y = d.model_KL;
bar(y);
hold on; p = scatter(1:length(y),y); set(p,'visible','off'); hold off
title(sprintf('Information gain (%d models)',nr),'FontSize',14);
xlabel('Dataset','FontSize',12); ylabel('Entropy (nats)','FontSize',12);

author: spm


text:
function [f,J,Q] = spm_fx_cmm(x,u,P,M)
% state equations for canonical neural-mass and mean-field models
% FORMAT [f,J,Q] = spm_fx_cmm(x,u,P,M)
%
% x - states and covariances
%
% x(i,j,k)        - k-th state of j-th population of i-th source
%                   i.e., running over sources, pop. and states
%
%   population: 1 - excitatory spiny stellate cells (input cells)
%               2 - superficial pyramidal cells     (forward output cells)
%               3 - inhibitory interneurons         (intrisic interneuons)
%               4 - deep pyramidal cells            (backward output cells)
%
%        state: 1 V  - voltage
%               2 gE - conductance (excitatory)
%               3 gI - conductance (inhibitory)
%
%--------------------------------------------------------------------------
% refs:
%
% Marreiros et al (2008) Population dynamics under the Laplace assumption
%
% See also:
%
% Friston KJ.
% The labile brain. I. Neuronal transients and nonlinear coupling. Philos
% Trans R Soc Lond B Biol Sci. 2000 Feb 29;355(1394):215-36.
%
% McCormick DA, Connors BW, Lighthall JW, Prince DA.
% Comparative electrophysiology of pyramidal and sparsely spiny stellate
% neurons of the neocortex. J Neurophysiol. 1985 Oct;54(4):782-806.
%
% Brunel N, Wang XJ.
% What determines the frequency of fast network oscillations with irregular
% neural discharges? I. Synaptic dynamics and excitation-inhibition
% balance. J Neurophysiol. 2003 Jul;90(1):415-30.
%
% Brunel N, Wang XJ.
% Effects of neuromodulation in a cortical network model of object working
% memory dominated by recurrent inhibition. J Comput Neurosci. 2001
% Jul-Aug;11(1):63-85.
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fx_cmm.m 5369 2013-03-28 20:09:27Z karl $

% get dimensions and configure state variables
%--------------------------------------------------------------------------
ns   = size(M.x,1);                      % number of sources
np   = size(M.x,2);                      % number of populations per source
nk   = size(M.x,3);                      % number of states per population
x    = reshape(x,ns,np,nk);              % hidden states

% extrinsic connection strengths
%==========================================================================

% exponential transform to ensure positivity constraints
%--------------------------------------------------------------------------
A{1} = exp(P.A{1});                      % forward
A{2} = exp(P.A{2});                      % backward
C    = exp(P.C);                         % subcortical

% detect and reduce the strength of reciprocal (lateral) connections
%--------------------------------------------------------------------------
for i = 1:length(A)
    L    = (A{i} > exp(-8)) & (A{i}' > exp(-8));
    A{i} = A{i}./(1 + 8*L);
end

% intrinsic connection strengths
%==========================================================================

% condition specific effects
%--------------------------------------------------------------------------
G    = full(P.H);
if any(P.G)
    G(2,2,:) = squeeze(G(2,2,:)) + P.G;
end
G(2,2,:) = 0;
G    = exp(G);

% connectivity switches
%==========================================================================
% 1 - excitatory spiny stellate cells (granular input cells)
% 2 - superficial pyramidal cells     (forward  output cells)
% 3 - inhibitory interneurons         (intrisic interneuons)
% 4 - deep pyramidal cells            (backward output cells)

% extrinsic connections (F B) - from superficial and deep pyramidal cells
%--------------------------------------------------------------------------
SA   = [1   0 ;
        0   1 ;
        0   2 ;
        0   0]/8;

% intrinsic connections (np x np) - excitatory
%--------------------------------------------------------------------------
GE   = [ 0    0    0    0
         4    0    0    0
         4    0    0    2
         0    4    0    0];

% intrinsic connections (np x np) - inhibitory
%--------------------------------------------------------------------------
GI   = [ 8    2    2    0
         0    8    2    0
         0    0   32    0
         0    0    8  128];

% rate constants (ns x np) (excitatory 4ms, inhibitory 16ms)
%--------------------------------------------------------------------------
KE   = 1000/4;                             % excitatory rate constants
KI   = exp(-P.T)*[1/64 1/32 1/16 1/16]*1000; % inhibitory rate constants

% Voltages
%--------------------------------------------------------------------------
VL   = -70;                               % reversal  potential leak (K)
VE   =  60;                               % reversal  potential excite (Na)
VI   = -90;                               % reversal  potential inhib (Cl)
VR   = -40;                               % threshold potential

CV   = exp(P.CV).*[16 16 32 16]/1000;     % membrane capacitance
GL   = 1;                                 % leak conductance

% mean-field effects:
%==========================================================================

% neural-mass approximation to covariance of states
%----------------------------------------------------------------------
Vx   = exp(P.S)*32;

% mean population firing and afferent extrinsic input
%--------------------------------------------------------------------------
m      = spm_Ncdf_jdw(x(:,:,1),VR,Vx);    % mean firing rate
a(:,1) = A{1}*m(:,2);                     % forward afference
a(:,2) = A{2}*m(:,4);                     % backward afference

% Averge background activity and exogenous input
%==========================================================================
BE     = exp(P.E)/2;

% input
%--------------------------------------------------------------------------
if isfield(M,'u')

    % endogenous input
    %----------------------------------------------------------------------
    U = u(:)*8/1000;

else

    % exogenous input
    %----------------------------------------------------------------------
    U = C*u(:)/1000;

end

% flow over every (ns x np) subpopulation
%==========================================================================
f     = x;
for i = 1:ns

        % intrinsic coupling
        %------------------------------------------------------------------
        E = (G(:,:,i).*GE)*m(i,:)';
        I = (G(:,:,i).*GI)*m(i,:)';

        % extrinsic coupling (excitatory only) and background activity
        %------------------------------------------------------------------
        E = E + BE + SA*a(i,:)';

        % and exogenous input(U)
        %------------------------------------------------------------------
        E(1) = E(1) + U(i);

        % Voltage
        %==================================================================
        f(i,:,1) =          (GL*(VL - x(i,:,1)) + ...
                      x(i,:,2).*(VE - x(i,:,1)) + ...
                      x(i,:,3).*(VI - x(i,:,1)) )./CV;

        % Conductance
        %==================================================================
        f(i,:,2) = (E' - x(i,:,2)).*KE;
        f(i,:,3) = (I' - x(i,:,3)).*KI(i,:);

end

% vectorise equations of motion
%==========================================================================
f = spm_vec(f);

if nargout < 2, return, end

% Jacobian
%==========================================================================
J = spm_cat(spm_diff(M.f,x,u,P,M,1));

if nargout < 3, return, end

% Delays
%==========================================================================
% Delay differential equations can be integrated efficiently (but
% approximately) by absorbing the delay operator into the Jacobian
%
%    dx(t)/dt     = f(x(t - d))
%                 = Q(d)f(x(t))
%
%    J(d)         = Q(d)df/dx
%--------------------------------------------------------------------------
% [specified] fixed parameters
%--------------------------------------------------------------------------
D  = [2 16];

d  = -D.*exp(P.D)/1000;
Sp = kron(ones(nk,nk),kron( eye(np,np),eye(ns,ns)));  % states: same pop.
Ss = kron(ones(nk,nk),kron(ones(np,np),eye(ns,ns)));  % states: same source

Dp = ~Ss;                            % states: different sources
Ds = ~Sp & Ss;                       % states: same source different pop.
D  = d(2)*Dp + d(1)*Ds;

% Implement: dx(t)/dt = f(x(t - d)) = inv(1 - D.*dfdx)*f(x(t))
%                     = Q*f = Q*J*x(t)
%--------------------------------------------------------------------------
Q  = spm_inv(speye(length(J)) - D.*J);


text:
function s = dicomLoadAllSeries(dicomDir, studyId, sortByFilenameFlag, ignoreAcquisition)
%
% s = dicomLoadAllSeries(dicomDir, [studyId=''], [sortByFilenameFlag=false])
%
% Loads all the dicom files found in the specified directory
% (recursively searches all sub-directories too). A structure
% array is returned, with one entry for each series found in the
% directory tree.
%
% If studyID (a string) is provided, then only series with a matching
% StudyID tag will be processed. Otherwise, all studies are processed.
%
% To save each series in a directory of DICOMs to separate NIFTI files, try
% niftiFromDicom.
%
% HISTORY:
% 2007.01.?? RFD: wrote it.
% 2009.12.08 RFD: completely replaced code to compute the imToScanner
% xform. I now use code that is modeled on NYU's dinifti, so we should now
% produce results consistent with that program. Which, BTW, are more
% accurate than what we did before.

if(~exist('sortByFilenameFlag','var')||isempty(sortByFilenameFlag))
    sortByFilenameFlag = false;
end
if(~exist('ignoreAcquisition','var')||isempty(ignoreAcquisition))
    ignoreAcquisition = false;
end

n = 0;
% Get a recursive dump of all sub-directories. We'll
% process every file in the dicomDir tree.
d = genpath(dicomDir);
if(isempty(d))
  error(['Dicom dir "' d '" not found or empty.']);
end
if(isunix)
    d = explode(':',genpath(dicomDir));
else
    d = explode(';',genpath(dicomDir));
end
d = d(1:end-1);
for(ii=1:length(d))
  d2 = dir(fullfile(d{ii},'*'));
  for(jj=1:length(d2))
	if(~d2(jj).isdir)
	  n = n+1;
	  dicomFiles{n} = fullfile(d{ii},d2(jj).name);
	end
  end
end

clear d;
numSeries = 0;
% We might need this if the dicoms are gzipped.
tmp = tempdir;
tmpName = tempname;

if(~exist('studyId','var'))
    studyId = [];
end
fprintf('Processing %d files in %s:\n', n, dicomDir);
for(imNum=1:length(dicomFiles))
    curFile = dicomFiles{imNum};
    if(length(curFile)>3 && strcmpi(curFile(end-2:end),'.gz'))
        if(exist(tmpName,'file')) delete(tmpName); end
        curFile = gunzip(curFile,tmp);
        movefile(curFile{1}, tmpName);
        curFile = tmpName;
    end
	try
	  info = dicominfo(curFile);
	  if((isempty(studyId) || strcmp(studyId,info.StudyID)) &&  isfield(info,'PatientPosition'))
        if(numSeries==0)
            curSeries = [];
        else
            if(ignoreAcquisition)
                curSeries = find(str2double(info.SeriesTime)==[s(:).acqTime]);
            else
                curSeries = find(info.AcquisitionNumber==[s(:).acqNum] & str2double(info.SeriesTime)==[s(:).acqTime]);
            end
        end
        if(isempty(curSeries))
		  numSeries = numSeries+1;
          curSeries = numSeries;
		  numSlice = 1;
		  s(curSeries).studyID = info.StudyID;
		  s(curSeries).studyDescription = getFieldVal(info, 'StudyDescription', '');
		  s(curSeries).studyDateTime = [info.SeriesDate ' ' info.StudyTime];
		  s(curSeries).patientName = info.PatientName.FamilyName;
		  s(curSeries).patientPosition = getFieldVal(info, 'PatientPosition', 'NONE');
		  s(curSeries).fieldStrength = getFieldVal(info, 'MagneticFieldStrength', []);
		  s(curSeries).seriesDescription = getFieldVal(info, 'SeriesDescription', []);
		  s(curSeries).acqTime = str2double(getFieldVal(info, 'SeriesTime', []));
		  s(curSeries).acqMatrix = getFieldVal(info, 'AcquisitionMatrix',[]);
		  s(curSeries).percentFOV = [getFieldVal(info, 'PercentSampling', []) getFieldVal(info, 'PercentPhaseFieldOfView', [])];
		  s(curSeries).reconDiam = getFieldVal(info, 'ReconstructionDiameter', []);
		  s(curSeries).sliceThickness = info.SliceThickness;
          if(isfield(info,'SpacingBetweenSlices'))
              s(curSeries).mmPerVox = [info.PixelSpacing(:)' info.SpacingBetweenSlices];
          else
              s(curSeries).mmPerVox = [info.PixelSpacing(:)' s(curSeries).sliceThickness];
          end
		  s(curSeries).TR = getFieldVal(info, 'RepetitionTime', []);
		  s(curSeries).TE = getFieldVal(info, 'EchoTime', []);
		  s(curSeries).inversionTime = getFieldVal(info, 'InversionTime', []);
		  s(curSeries).SAR = getFieldVal(info, 'SAR', []);
		  s(curSeries).pixBandwidth = getFieldVal(info, 'PixelBandwidth', []);
          if(~isempty(s(curSeries).pixBandwidth) && numel(s(curSeries).acqMatrix>=2))
              s(curSeries).Bandwidth = s(curSeries).pixBandwidth*s(curSeries).acqMatrix(2);
          else
              s(curSeries).Bandwidth = [];
          end
		  s(curSeries).NEX = getFieldVal(info, 'NumberOfAverages', []);
		  s(curSeries).imageFreq = getFieldVal(info, 'ImagingFrequency',[]);
		  s(curSeries).flipAngle = getFieldVal(info, 'FlipAngle', []);
		  s(curSeries).phaseEncodeDir = getFieldVal(info, 'InPlanePhaseEncodingDirection', []);
		  % We'll fill in the 3rd and 4th dims below
		  %s(curSeries).dims = [info.Rows info.Columns info.ImagesInAcquisition];
		  s(curSeries).dims = [info.Rows info.Columns 0 0];
		  s(curSeries).seriesNum = getFieldVal(info, 'SeriesNumber', []);
		  s(curSeries).acqNum = getFieldVal(info, 'AcquisitionNumber', []);
		  s(curSeries).imageOrientation = getFieldVal(info, 'ImageOrientationPatient', []);
          if(strcmp(info.Manufacturer,'GE MEDICAL SYSTEMS'))
              % *** TO DO: these should be pulse-sequence dependent
              s(curSeries).sequenceName = getFieldVal(info, 'Private_0019_109e',[]);
              s(curSeries).mtOffset = getFieldVal(info, 'Private_0043_1034',[]);
              s(curSeries).dtiBValue = str2double(char(getFieldVal(info, 'Private_0019_10b0',[])));
              s(curSeries).dtiGradsCode = str2double(char(getFieldVal(info, 'Private_0019_10b2',[])));
          elseif(strcmp(info.Manufacturer,'SIEMENS'))
              bmtx = getFieldVal(info, 'Private_0019_1027', []);
              if(isempty(bmtx))
                  s(curSeries).dtiBMatrix = [];
                  s(curSeries).dtiBValue = 0;
                  s(curSeries).dtiGradDir = [1 0 0]';
              else
                  s(curSeries).dtiBMatrix = typecast(bmtx,'double');
                  s(curSeries).dtiBValue = norm(sqrt(s(curSeries).dtiBMatrix([1 4 6])))^2;
                  s(curSeries).dtiGradDir = sign(sign(s(curSeries).dtiBMatrix([1:3])) + 0.01).*sqrt(s(curSeries).dtiBMatrix([1 4 6])/s(curSeries).dtiBValue);
              end
          end
		  %ReferencedImageSequence

		  %s(curSeries).imData = zeros(s(curSeries).dims);
		  %s(curSeries).imagePosition = zeros(3,s(curSeries).dims(3));
		  %s(curSeries).sliceNum = zeros(1,s(curSeries).dims(3));
		  % We'll set this to more meaningful values below
		  s(curSeries).imToScanXform = eye(4);

		  fprintf('Loading series %d, acquisition %d...\n', s(curSeries).seriesNum, s(curSeries).acqNum);
        end
        s(curSeries).sliceNum(numSlice) = info.InstanceNumber;
        s(curSeries).imagePosition(:,numSlice) = info.ImagePositionPatient;
        s(curSeries).sliceLoc(numSlice) = info.SliceLocation;
        s(curSeries).fileName{numSlice} = dicomFiles{imNum};
        % Matlab seems to flip permute x and y when reading the data, so we do
        % a transpose here. Seems to make everything work better down the line.
        % *** RFD: need to flip dim 2 also?
        s(curSeries).imData(:,:,numSlice) = dicomread(info)';
        numSlice = numSlice+1;
	  end
	catch
	  disp(['failed to load ' curFile]);
	end
end

% sort the series properly
[junk,ind] = sort([s(:).seriesNum]);
s = s(ind);

% Sort the slices based on sliceNum and break timeseries up into separate
% volumes (good for fMRI and DTI acquisitions).

for(ii=1:numSeries)
    if(sortByFilenameFlag)
        [s(ii).sliceNum,sortInd] = sort(s(ii).fileName);
    else
        % Sort based on sliceNum (DICOM InstanceNumber).
        [s(ii).sliceNum,sortInd] = sort(s(ii).sliceNum);
    end
    s(ii).sliceLoc = s(ii).sliceLoc(sortInd);
    s(ii).imagePosition = s(ii).imagePosition(:,sortInd);
    s(ii).imData = s(ii).imData(:,:,sortInd);
    % Now compute the actual number of slices per timepoint
    uniqueSliceLoc = unique(s(ii).sliceLoc);
    s(ii).dims(3) = length(uniqueSliceLoc);
    s(ii).dims(4) = floor(length(sortInd)./s(ii).dims(3));
    if(prod(double(s(ii).dims))==numel(s(ii).imData))
        s(ii).imData = reshape(s(ii).imData, s(ii).dims);
    else
        warning(['Series ' num2str(ii) ' appears to contain >1 volume, but I can not guess the dimensions- keeping it a single vol.']);
    end
end

for(ii=1:length(s))
    if(isempty(s(ii).imData)) s(ii) = []; end
end
numSeries = length(s);

% Compute the scanner-to-image xform and resort the slices based on the
% slice normal.
for(ii=1:numSeries)
    % 2009.12.07 RFD: Pulled orientation code from NYU's dinifti
    dim = size(s(ii).imData);

    % Orientation information
    %-------------------------------------------------------------------
    % Axial Analyze voxel co-ordinate system:
    % x increases     right to left
    % y increases posterior to anterior
    % z increases  inferior to superior

    % DICOM patient co-ordinate system:
    % x increases     right to left
    % y increases  anterior to posterior
    % z increases  inferior to superior

    % T&T co-ordinate system:
    % x increases      left to right
    % y increases posterior to anterior
    % z increases  inferior to superior

    % Code lifted from dinifti source:
    % patient orientation cosines, in [Sag, Cor, Ax] order
    rowCosines = s(ii).imageOrientation(1:3);
    colCosines = s(ii).imageOrientation(4:6);
    sliceNorm = cross(rowCosines, colCosines);

    qto_xyz = zeros(4,4);
    qto_xyz(1,1) = -rowCosines(1); % -image->RowSagCos();
	qto_xyz(1,2) = -colCosines(1); % -image->ColSagCos();
	qto_xyz(1,3) = -sliceNorm(1); % -(rowCosines(2) * colCosines(3) - rowCosines(3) * colCosines(2)); % -image->SagNorm();

	qto_xyz(2,1) = -rowCosines(2); % -image->RowCorCos();
	qto_xyz(2,2) = -colCosines(2); % -image->ColCorCos();
	qto_xyz(2,3) = -sliceNorm(2); % -(rowCosines(3) * colCosines(1) - rowCosines(1) * colCosines(3)); % -image->CorNorm();

	qto_xyz(3,1) = rowCosines(3); % image->RowTraCos();
	qto_xyz(3,2) = colCosines(3); % image->ColTraCos();
	qto_xyz(3,3) = sliceNorm(3); % rowCosines(1) * colCosines(2) - rowCosines(2) * colCosines(1); % image->TraNorm();

    % *** TO DO: The most robust wayt to define slice order and slice
    % spacing is to project slice position onto slice normal.
    % From: http://nifti.nimh.nih.gov/board/read.php?f=1&i=579&t=508
    % Slice normal vector is the cross product of ImageOrientationPatient vectors.
    % Project each slice's ImagePositionPatient projection onto the slice normal
    % to find it's position in the slice order list.
    % pfirst = dot(sliceNorm,s(ii).imagePosition(:,1)) * sliceNorm;
    % plast = dot(sliceNorm,s(ii).imagePosition(:,end)) * sliceNorm;

    % Simple hack to decide if we need to flip the slice order:
    if(dot(sliceNorm,s(ii).imagePosition(:,1)) > dot(sliceNorm,s(ii).imagePosition(:,end)))
        s(ii).sliceNum = flipdim(s(ii).sliceNum,2);
        s(ii).sliceLoc = flipdim(s(ii).sliceLoc,2);
        s(ii).imagePosition = flipdim(s(ii).imagePosition,2);
        s(ii).imData = flipdim(s(ii).imData,3);
    end

	% offset = slicePosition(center) - M*(i,j,k of center) since we are working with the first slice k = 0.

    % For matlab's 1-indexing, we have to subtract mmPerVox from the origin
    % (imagePosition) to reflect the fact that the first voxel is [1,1,1]
    % rather than [0,0,0].
    pos = s(ii).imagePosition(:,1);
    qto_xyz(:,4) = [-pos(1) -pos(2) pos(3) 1]';
    qto_xyz(1:3,1:3) = qto_xyz(1:3,1:3)*diag(s(ii).mmPerVox);

    % Apply a 1-voxel offset to the origin to account for Matlab's
    % 1-indexing.
    qto_xyz = inv(inv(qto_xyz)+[0 0 0 1;0 0 0 1;0 0 0 1; 0 0 0 0]);

    s(ii).imToScanXform  = qto_xyz;

% OLD CODE:
%     s(ii).imToScanXform = eye(4);
%     % DICOM +Xd is Left, +Yd is Posterior, +Zd is Superior,
%     % while NIFTI is +x is Right, +y is Anterior, +z is Superior.
%     % So, x and y offsets get flipped.
%     %s(ii).imToScanXform(1:3,4) = s(ii).imagePosition(:,1).*[-1 -1 1]';
%     s(ii).imToScanXform(1:3,4) = -s(ii).imagePosition(:,1);
%
%     % From the NIFTI-1 standard (Bob Cox):
%     % The DICOM attribute (0020,0037) "Image Orientation (Patient)"
%     % gives the orientation of the x- and y-axes of the image data
%     % in terms of 2 3-vectors. The first vector is a unit vector
%     % along the x-axis, and the second is along the y-axis. If the
%     % (0020,0037) attribute is extracted into the value
%     % (xa,xb,xc,ya,yb,yc), then the first two columns of the R
%     % matrix would be
%     %            [ -xa  -ya ]
%     %            [ -xb  -yb ]
%     %            [  xc   yc ]
%     % The negations are because DICOM's x- and y-axes are reversed
%     % relative to NIFTI's. The third column of the R matrix gives
%     % the direction of displacement (relative to the subject) along
%     % the slice-wise direction. This orientation is not encoded in
%     % the DICOM standard in a simple way; DICOM is mostly concerned
%     % with 2D images.  The third column of R will be either the
%     % cross-product of the first 2 columns or its negative.  It is
%     % possible to infer the sign of the 3rd column by examining the
%     % coordinates in DICOM attribute (0020,0032) "Image Position
%     % (Patient)" for successive slices.  However, this method
%     % occasionally fails for reasons that I (RW Cox) do not understand.
%     s(ii).imToScanXform(1:3,1) = s(ii).imageOrientation(1:3).*[-1 -1 1]';
%     s(ii).imToScanXform(1:3,2) = s(ii).imageOrientation(4:6).*[-1 -1 1]';
%     s(ii).imToScanXform(1:3,3) = cross(s(ii).imToScanXform(1:3,1),s(ii).imToScanXform(1:3,2));
%     if(s(ii).sliceLoc(1)>s(ii).sliceLoc(end))
%         % Slices run S-I, so signal a flip in this dim.
%         s(ii).imToScanXform(1:3,3) = -s(ii).imToScanXform(1:3,3);
%         s(ii).imToScanXform(3,4) = -s(ii).imToScanXform(3,4);
%     end
%
%     % CHECK THIS: How to infer the orientation?
%     % Check to see if the image plane normal goes away from the
%     % center of the volume?
%     %if(s(ii).slicePos(1)<s(ii).slicePos(end))
%     %  s(ii).imToScanXform(1:3,3) = -s(ii).imToScanXform(1:3,3);
%     %end
%     s(ii).imToScanXform = s(ii).imToScanXform*diag([s(ii).mmPerVox 1]);
end

return;

function val = getFieldVal(s, fieldName, defaultValue)
    if(isfield(s,fieldName))
        val = s.(fieldName);
    else
        val = defaultValue;
    end
return;


text:
function fx = camel ( n, x )

%*****************************************************************************80
%
%% CAMEL evaluates the six-hump camel back function.
%
%  Discussion:
%
%    This function has been modified from its original definition by
%    multiplying by -1 and then adding 1, so that there is a small island
%    of positivity within the box [-2.5,2.5]x[-2.5,2.5].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N,2), the point coordinates.
%
%    Output, real FX(N), the function values.
%
  fx = zeros ( n, 1 );

  fx = - ...
  ( ...
      ( 4 - 2.1 * x(:,1).^2 + 0.33 * x(:,1).^4 ) .* x(:,1).^2 ...
    + x(:,1) .* x(:,2) + 4 * ( x(:,2).^2 - 1 ) .* x(:,2).^2 ...
  ) + ones ( n, 1 );

  return
end


text:
% function [Tk,dTk] = interp2DlinearP(TD,Omega,X,doderivative)
%
% (c) NP, Safir, Luebeck 2006
% this is a function for linear interpolation
% input:
%     - TD      :   the data that has to be interpolated (a 3D matlab
%     matrix)
%     - Omega   :   the image domain
%     - X       :   the interpolaiton points (a long vector; see getGrid
%     for example)
%     - doderivative : flag for computing the derivative
%
% output:
%     - Tk      : the interpolated function values (same format as X)
%     - dTk     : the derivative of the linear interpolated function with
%     respect to X (computed when doderivative is 1)
%
function [Tk,dTk] = interp2DlinearP(TD,Omega,X,periode,doderivative);

Tk = [];
dTk = [];

% if no input is given a testcase will be created and then
% the derivative test will be started
if nargin == 0
  fprintf('%s\n',mfilename);
  Omega = [1, 1];
  %TD = [zeros(2,7); ones(3,1)*[0 0 0.5 1 1.5 0 0.5];zeros(2,7)];
  TD = [0 1.5;  0 1.5];
  TD = TD';
  periode = 2;
%   figure(1);clf;
%   image(TD./periode.*256); colormap(hsv(256)); axis image;

  m = [100,100];
  X = getGrid(Omega,m);
  %    m = [1,1];
  %    X = [0.5; 0.5];
interp2DlinearP(TD,Omega,X,periode,0);

%    viewImage(Tk./periode*256,Omega,m,'fig',2);
%     figure(3); clf;
%     surf(reshape(Tk,m)');
  return;
end;
% if nargout == 0 start a derivative test
if nargout == 0
  % derivative test
  fprintf('derivative(%s)\n',mfilename);
  [fc,dfc] = feval(mfilename,TD,Omega,X,periode,1);
  v   = randn(size(X));
  %   %save v1 v
  %   load v1
  %   e = zeros(size(v));
  %   e(1) = v(1);
  %   v = e;
  %   v
  dfv = dfc*v;
  hh  = logspace(0,-10,11);
  fprintf('%12s %12s %12s \n','h','|fc-ft|','|fc+vdfc-ft|');
  for j=1:length(hh),
    Xt = X + hh(j)*v;
    ft = feval(mfilename,TD,Omega,Xt,periode,0);
    n1 = norm(fc(:)-ft(:));
    n2 = norm(fc(:)+hh(j)*dfv-ft(:));
    fprintf('%12.4e %12.4e %12.4e\n',hh(j),n1,n2);
  end;
  return;
end;

% if doderivative is not explicitly given
% a derivative is computed if two output arguments are given
if isempty(doderivative), doderivative = (nargout > 0); end;
% n is the number of interpolation points
% mD is the size of the interpolation dat
n = length(X)/2;
mD = size(TD);
periode2 = periode/2;
% zeropadding to avoid bad cases
% so we do an embedding of TD into a field of zeros
offset = 1;
To = zeros(size(TD)+2*offset);
To(offset+(1:size(TD,1)),offset+(1:size(TD,2))) = TD;

% get pixelsize, pay attention to the change of order from m to mD
hD = Omega./mD([2,1]);
% for easier reading
hD1 = hD(1);
hD2 = hD(2);
% for easier reading

% alloc memory for output

% alloc memory for output
Tk = zeros(n,1);
if doderivative, dTk = zeros(n,2); end;

X1 = X(1:n); X2 = X(n+(1:n));

% transform grid X to integer grid
X1 = (1/hD1)*X1 + 0.5;
X2 = (1/hD2)*X2 + 0.5;

% so if we have no good points, we have to format the
% derivative for numerical reasons

G = find( 0 < X1 & X1 < mD(2) + 1 & ...
  0 < X2 & X2 < mD(1) + 1    );

if isempty(G),
  dTk = sparse(n,2*n);
  return;
end;

% dj is the size we have to go to get the
% neigbhour in j direction
d1 = size(To,1);
d2 = -1;

% now we get into more detail
% Kj storages the noninteger part of Xj
K1 = floor(X1(G)); K2 = floor(X2(G));
Xi1 = X1(G) - K1;  Xi2 = X2(G) - K2;

K = offset + (offset + K1 - 1)*d1 + (mD(1) - K2 +1);

% compute Tk
% here now something new
dTs1 = To(K+d1)    - To(K);
dTs2 = To(K+d1+d2) - To(K+d2);

pO1 = zeros(size(K));
pO1(dTs1 > periode2)  = -periode;
pO1(dTs1 < -periode2) =  periode;

pO2 = zeros(size(K));
pO2(dTs2 > periode2)  = -periode;
pO2(dTs2 < -periode2) =  periode;

Tk1 = zeros(size(G));
Tk2 = zeros(size(G));

Tk1 =  To(K) .* (1-Xi1)  +  (To(K+d1) + pO1) .* Xi1;
Tk1(Tk1<0)        = Tk1(Tk1<0)        + periode;
Tk1(Tk1>=periode) = Tk1(Tk1>=periode) - periode;

Tk2 =  To(K+d2) .* (1-Xi1)  +  (To(K+d1+d2) + pO2) .* Xi1;
Tk2(Tk2<0)        = Tk2(Tk2<0)        + periode;
Tk2(Tk2>=periode) = Tk2(Tk2>=periode) - periode;

dTs3 = Tk2 - Tk1;
pO3 = zeros(size(K));
pO3(dTs3 > periode2)  = -periode;
pO3(dTs3 < -periode2) =  periode;

%keyboard;
Tk(G) = Tk1 .* (1-Xi2)  +  (Tk2 + pO3) .* Xi2;

Tk(Tk<0)        = Tk(Tk<0)        + periode;
Tk(Tk>=periode) = Tk(Tk>=periode) - periode;

% keyboard

% if derivative is needed compute it
% next step: the derivative ???

if doderivative,
  dTk(G,1) = (To(K+d1) + pO1 - To(K)   ) .* (1-Xi2) ...
    + (To(K+d1+d2) + pO2 - To(K+d2)).* (Xi2);
  dTk(G,2) = Tk2 + pO3 - Tk1;
end;
if doderivative,
  dTk = spdiags([dTk(:,1)/hD1, dTk(:,2)/hD2],[0,n],n,2*n);
end;

return;


text:
function [varargout] = likUni(hyp, y, mu, s2, inf, i)

% likUni - Uniform likelihood function for classification. The expression for
% the likelihood is
%   likUni(t) = 1/2.
%
% There are no hyperparameters:
%
% hyp = [ ]
%
% Several modes are provided, for computing likelihoods, derivatives and moments
% respectively, see likFunctions.m for the details. In general, care is taken
% to avoid numerical issues when the arguments are extreme.
%
% Copyright (c) by Hannes Nickisch, 2013-09-02.
%
% See also LIKFUNCTIONS.M.

if nargin<3, varargout = {'0'}; return; end   % report number of hyperparameters

if nargin<5                              % prediction mode if inf is not present
  lp = -log(2)*ones(size(mu));
  ymu = {}; ys2 = {};
  if nargout>1
    p = exp(lp);
    ymu = 2*p-1;                                                % first y moment
    if nargout>2
      ys2 = 4*p.*(1-p);                                        % second y moment
    end
  end
  varargout = {lp,ymu,ys2};
else
  switch inf
  case 'infLaplace'
    if nargin<6                                             % no derivative mode
      lp = -log(2)*ones(size(mu)); dlp = {}; d2lp = {}; d3lp = {};
      if nargout>1
        dlp = zeros(size(mu));               % dlp, derivative of log likelihood
        if nargout>2                    % d2lp, 2nd derivative of log likelihood
          d2lp = zeros(size(mu));
          if nargout>3                  % d3lp, 3rd derivative of log likelihood
            d3lp = zeros(size(mu));
          end
        end
      end
      varargout = {lp,dlp,d2lp,d3lp};
    else                                                       % derivative mode
      varargout = {[],[],[]};                         % derivative w.r.t. hypers
    end

  case 'infEP'
    if nargin<6                                             % no derivative mode
      lZ = -log(2)*ones(size(mu));                           % log part function
      dlZ = {}; d2lZ = {};
      if nargout>1
        dlZ  = zeros(size(mu));                     % 1st derivative w.r.t. mean
        if nargout>2
          d2lZ = zeros(size(mu));                   % 2nd derivative w.r.t. mean
        end
      end
      varargout = {lZ,dlZ,d2lZ};
    else                                                       % derivative mode
      varargout = {[]};                                     % deriv. wrt hyp.lik
    end

  case 'infVB'
    % variational lower site bound
    % t(s) = 1/2
    % the bound has the form: (b+z/ga)*f - f.^2/(2*ga) - h(ga)/2
    n = numel(s2); b = zeros(n,1); y = y.*ones(n,1); z = y;
    varargout = {b,z};
  end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [pdag, G] = learn_struct_pdag_pc(cond_indep, n, k, varargin)
% LEARN_STRUCT_PDAG_PC Learn a partially oriented DAG (pattern) using the PC algorithm
% P = learn_struct_pdag_pc(cond_indep, n, k, ...)
%
% n is the number of nodes.
% k is an optional upper bound on the fan-in (default: n)
% cond_indep is a boolean function that will be called as follows:
%   feval(cond_indep, x, y, S, ...)
% where x and y are nodes, and S is a set of nodes (positive integers),
% and ... are any optional parameters passed to this function.
%
% The output P is an adjacency matrix, in which
% P(i,j) = -1 if there is an i->j edge.
% P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j
%
% The PC algorithm does structure learning assuming all variables are observed.
% See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117.
% This algorithm may take O(n^k) time if there are n variables and k is the max fan-in,
% but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n).

sep = cell(n,n);
ord = 0;
done = 0;
G = ones(n,n);
G=setdiag(G,0);
while ~done
  done = 1;
  [X,Y] = find(G);
  for i=1:length(X)
    x = X(i); y = Y(i);
    %nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);
    nbrs = mysetdiff(neighbors(G, y), x);  % bug fix by Raanan Yehezkel <raanany@ee.bgu.ac.il> 6/27/04
    nbrs = unique(nbrs); % bug fix Immer Ebert ebert@tree.com 3/278/11
    if length(nbrs) >= ord & G(x,y) ~= 0
      done = 0;
      %SS = subsets(nbrs, ord, ord); % all subsets of size ord
      SS = subsets1(nbrs, ord);
      for si=1:length(SS)
	S = SS{si};
	if feval(cond_indep, x, y, S, varargin{:})
	  %if isempty(S)
 	  %  fprintf('%d indep of %d ', x, y);
	  %else
 	  %  fprintf('%d indep of %d given ', x, y); fprintf('%d ', S);
 	  %end
 	  %fprintf('\n');

	  % diagnostic
	  %[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:});
	  %fprintf(': r = %6.4f\n', r);

	  G(x,y) = 0;
	  G(y,x) = 0;
	  sep{x,y} = myunion(sep{x,y}, S);
	  sep{y,x} = myunion(sep{y,x}, S);
	  break; % no need to check any more subsets
	end
      end
    end
  end
  ord = ord + 1;
end

% Create the minimal pattern,
% i.e., the only directed edges are V structures.
pdag = G;
[X, Y] = find(G);
% We want to generate all unique triples x,y,z
% This code generates x,y,z and z,y,x.
for i=1:length(X)
  x = X(i);
  y = Y(i);
  Z = find(G(y,:));
  Z = mysetdiff(Z, x);
  for z=Z(:)'
    if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x})
      %fprintf('%d -> %d <- %d\n', x, y, z);
      pdag(x,y) = -1; pdag(y,x) = 0;
      pdag(z,y) = -1; pdag(y,z) = 0;
    end
  end
end

% Convert the minimal pattern to a complete one,
% i.e., every directed edge in P is compelled
% (must be directed in all Markov equivalent models),
% and every undirected edge in P is reversible.
% We use the rules of Pearl (2000) p51 (derived in Meek (1995))

old_pdag = zeros(n);
iter = 0;
while ~isequal(pdag, old_pdag)
  iter = iter + 1;
  old_pdag = pdag;
  % rule 1
  [A,B] = find(pdag==-1); % a -> b
  for i=1:length(A)
    a = A(i); b = B(i);
    C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a
    if ~isempty(C)
      pdag(b,C) = -1; pdag(C,b) = 0;
      %fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C);
    end
  end
  % rule 2
  [A,B] = find(pdag==1); % unoriented a-b edge
  for i=1:length(A)
    a = A(i); b = B(i);
    if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' );
      pdag(a,b) = -1; pdag(b,a) = 0;
      %fprintf('rule 2: %d -> %d\n', a, b);
    end
  end
  % rule 3
  [A,B] = find(pdag==1); % a-b
  for i=1:length(A)
    a = A(i); b = B(i);
    C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' );
    % C contains nodes c s.t. a-c->ba
    G2 = setdiag(G(C, C), 1);
    if any(G2(:)==0) % there are 2 different non adjacent elements of C
      pdag(a,b) = -1; pdag(b,a) = 0;
      %fprintf('rule 3: %d -> %d\n', a, b);
    end
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [clusteredEdges, clusterInfo] = clusterCandidatePairs_3logFeatures_normalized_kmeans(imdb, imageSet, varargin)
%clusterCandidatePairs_3logFeatures_normalized_kmeans assigns all pairs of candidates to clusters, k-means clustering is performed if necessary

if ~exist('varargin', 'var')
    varargin = {};
end
%% parameters
opts = struct;
opts.numClusters = 10;
opts.clusterInfo = struct('type', [], 'mean', [], 'std', [], 'numClusters', [], 'clusterCenters', [] );
opts.maxPointsToCluster = 10^6;
opts.labelIouThreshold = 0.4;
opts.scoreThreshold = 0;
opts.dataPath = '';
opts.candidateIds = [];
opts.randomSeed = 1;
opts = vl_argparse(opts, varargin);

if ~exist('imageSet', 'var') || isempty( imageSet )
    error( 'clusterCandidatePairs_3logFeatures_normalized_kmeans:emptyImageSet', 'There are not images to do clustering' );
end

%% get the candidate data
dataset = struct;
dataset.boundingBoxes = cell(length(imageSet), 1);
dataset.numCandidates = zeros(length(imageSet), 1);
dataset.numCandidatesInitial = zeros(length(imageSet), 1);
dataset.groundTruth = cell(length(imageSet), 1);
dataset.scores = cell(length(imageSet), 1);
dataset.imageSize = cell(length(imageSet), 1);

numBoxPairs = 0;

selectedCandidatesIds = cell(length(imageSet), 1);
for iImageId = 1 : length(imageSet)
    if mod( iImageId, 1000 ) == 0
        fprintf('Reading data for image %d of %d\n', iImageId, length(imageSet));
    end

    iImage = imageSet( iImageId );

    curCandidates = load( fullfile( opts.dataPath, imdb.candidateFiles{ iImage } ) );

    dataset.numCandidatesInitial(iImageId) = size( curCandidates.boxes, 1);

    selectedCandidatesIds{iImageId} = [];
    if isempty( opts.candidateIds )
        selectedCandidatesIds{iImageId} = imdb.candidateIds{ iImage };
    else
        selectedCandidatesIds{iImageId} = opts.candidateIds{iImageId};
    end
    curCandidates = curCandidates.boxes( selectedCandidatesIds{iImageId}, :);

    % fix the Bb format: SelectiveSearch format [y1 x1 y2 x2] to format [x y w h]
    curBb = convertBb_Y1X1Y2X2_to_X1Y1WH( curCandidates(:, 1 : 4) );

    % fix the ground-truth format
    curGt = convertBb_X1Y1X2Y2_to_X1Y1WH( imdb.groundTruth{iImage} );

    dataset.imageSize{iImageId} = imdb.images.size{iImage};
    dataset.numCandidates(iImageId) = size(curBb, 1);
    dataset.boundingBoxes{iImageId} = curBb;
    dataset.groundTruth{iImageId} = curGt;
    dataset.scores{iImageId} = imdb.scores{iImage};

    numBoxPairs = numBoxPairs + size(curBb, 1) * (size(curBb, 1) - 1) / 2;
end

%% get info
numObjects = length(dataset.numCandidates);
numBoxFeatures = 8;
bBoxData = nan(numBoxPairs, numBoxFeatures);
bBoxPairLabels = nan(numBoxPairs, 2);
bBoxPairScores = nan(numBoxPairs, 2);
image_w = nan(numBoxPairs, 1);
image_h = nan(numBoxPairs, 1);
imageId = nan(numBoxPairs, 1);

clusteredEdges = cell(numObjects, 1);
startObjectId = zeros(numObjects + 1, 1);

%% collect all pairs
iPair = 0;
for iObject = 1 : numObjects
    startObjectId(iObject) = iPair + 1;
    curBbNum = dataset.numCandidates(iObject);
    curPairNum = curBbNum * (curBbNum - 1) / 2;
    curIds = iPair + 1 : iPair + curPairNum;
    iPair = iPair + curPairNum;

    [bb1, bb2] = meshgrid(1 : curBbNum, 1 : curBbNum);

    curMask = bb1 < bb2;
    bb1 = bb1(curMask);
    bb1 = bb1(:);
    bb2 = bb2(curMask);
    bb2= bb2(:);

    curBBoxes = dataset.boundingBoxes{ iObject };

    curGtIoU = zeros( size(curBBoxes, 1), 1 );
    for iGt = 1 : size( dataset.groundTruth{iObject}, 1 )
        IoU = bbIntersectionOverUnion( curBBoxes, dataset.groundTruth{iObject}(iGt, :) );
        curGtIoU = max( curGtIoU, IoU );
    end
    curLabels = curGtIoU > opts.labelIouThreshold;

    bBoxData(curIds, :) = [ curBBoxes( bb1, : ), curBBoxes( bb2, : ) ];
    bBoxPairLabels(curIds, :) = [ curLabels(bb1), curLabels(bb2) ];
    bBoxPairScores(curIds, :) = [ dataset.scores{iObject}(bb1), dataset.scores{iObject}(bb2) ];

    image_w( curIds ) = dataset.imageSize{iObject}( 2 );
    image_h( curIds ) = dataset.imageSize{iObject}( 1 );
    imageId( curIds ) = iObject;

    clusteredEdges{iObject} = struct;
    clusteredEdges{iObject}.bbIds = [bb1, bb2];
    clusteredEdges{iObject}.clusterId = nan(size(bb1, 1), 1);
end
startObjectId( numObjects + 1 ) = iPair + 1;
bBoxData(iPair + 1 : end, :) = [];
image_w(iPair + 1 : end, :) = [];
image_h(iPair + 1 : end, :) = [];
bBoxPairLabels(iPair + 1 : end, :) = [];
bBoxPairScores(iPair + 1 : end, :) = [];
imageId(iPair + 1 : end, :) = [];

%% sort bouding boxes (left to right)
changeOrder = ( bBoxData(:, 1) > bBoxData(:, 5) ) | ( (bBoxData(:, 1) == bBoxData(:, 5)) & (bBoxData(:, 2) > bBoxData(:, 6)) );
tmp = bBoxData(changeOrder, 1 : 4);
bBoxData(changeOrder, 1 : 4) = bBoxData(changeOrder, 5 : 8);
bBoxData(changeOrder, 5 : 8) = tmp;

tmp = bBoxPairLabels(changeOrder, 1);
bBoxPairLabels(changeOrder, 1) = bBoxPairLabels(changeOrder, 2);
bBoxPairLabels(changeOrder, 2) = tmp;

tmp = bBoxPairScores(changeOrder, 1);
bBoxPairScores(changeOrder, 1) = bBoxPairScores(changeOrder, 2);
bBoxPairScores(changeOrder, 2) = tmp;

%% get BB features
cx1 = bBoxData(:, 1) + bBoxData(:, 3) / 2;
cx2 = bBoxData(:, 5) + bBoxData(:, 7) / 2;
cy1 = bBoxData(:, 2) + bBoxData(:, 4) / 2;
cy2 = bBoxData(:, 6) + bBoxData(:, 8) / 2;
w1 = bBoxData(:, 3);
w2 = bBoxData(:, 7);
h1 = bBoxData(:, 4);
h2 = bBoxData(:, 8);
s1 = (w1 + h1) / 2;
s2 = (w2 + h2) / 2;

f1 = s1 ./ s2;
f2 = (cx2 - cx1) ./ s1;
f3 = (cy2 - cy1) ./ s1;
f1 = log(f1);
f2 = log( abs(f2) + 1 ) .* sign(f2);
f3 = log( abs(f3) + 1 ) .* sign(f3);
features = [f1, f2, f3];

%% clustering
if isempty( opts.clusterInfo.type )
    %% Run kmeans
    rng(opts.randomSeed , 'twister');

    % normalize the features
    meanFeatures = mean( features, 1 );
    stdFeatures = std( features, 0, 1 );

    featuresNormalized = bsxfun( @minus, features, meanFeatures );
    featuresNormalized = bsxfun( @rdivide, featuresNormalized, stdFeatures );

    clusterInfo = struct;
    clusterInfo.type = '3logFeatures_normalized';
    clusterInfo.mean = meanFeatures;
    clusterInfo.std = stdFeatures;
    clusterInfo.numClusters = opts.numClusters;
    clusterInfo.clusterCenters = [];

    curMask = ( ( bBoxPairScores(:, 1) > opts.scoreThreshold ) | ( bBoxPairLabels(:, 1) == 1) ) ...
        & ( ( bBoxPairScores(:, 2) > opts.scoreThreshold ) | ( bBoxPairLabels(:, 2) == 1) );

    curIds = find( curMask );
    if numel( curIds ) > opts.maxPointsToCluster
        curIds = curIds( randperm( numel( curIds ), opts.maxPointsToCluster ) );
    end

    [~,  clusterInfo.clusterCenters] = kmeans( featuresNormalized(curIds, :), opts.numClusters );

else
    clusterInfo = opts.clusterInfo;
    featuresNormalized = bsxfun( @minus, features, clusterInfo.mean );
    featuresNormalized = bsxfun( @rdivide, featuresNormalized, clusterInfo.std );
end
clusterIds = assignPointsToClusters( clusterInfo.clusterCenters, featuresNormalized );

%% change order within some pairs bounding boxes, revert to original indices
for iObject = 1 : numObjects
    curIds = startObjectId(iObject) : 1 : startObjectId(iObject + 1) - 1;
    changeFlags = changeOrder( curIds );
    curClusterIds = clusterIds( curIds );

    tmp = clusteredEdges{iObject}.bbIds(changeFlags, 1);
    clusteredEdges{iObject}.bbIds(changeFlags, 1) = clusteredEdges{iObject}.bbIds(changeFlags, 2);
    clusteredEdges{iObject}.bbIds(changeFlags, 2) = tmp;
    clusteredEdges{iObject}.clusterId = curClusterIds;

    clusteredEdges{iObject}.initialBbIds = selectedCandidatesIds{ iObject }( clusteredEdges{iObject}.bbIds );
end

end

author: aosokin


text:
function track(imagefig,varargins)
global M
global Nsamples

%%% Setting figure properties to start tracking
set(gcf,'windowbuttondownfcn',@dataout);
set(gcf,'WindowButtonMotionFcn',@datain);
set(gcf,'userdata',[]);

%%% Creating initial input, since the filter will need samplea at
%%% n=0,-1,-2... When processing the X[n=1] sample.
for k=1:Nsamples
    temp=get(gca,'currentpoint');
    set(gcf,'userdata',[get(gcf,'userdata'); [temp(1,1:2), 0 ]]);
end
tic; % Begin time measuring


text:
function tests = test_spx_norms
  tests = functiontests(localfunctions);
end

function test_lp_norms(testCase)
    x = [1 1 1
         2 2 2
        3 3 3];
    verifyEqual(testCase, [6 6 6], spx.norm.norms_l1_cw(x));
    verifyEqual(testCase, [3 6 9]', spx.norm.norms_l1_rw(x));
    verifyEqual(testCase, sqrt([14 14 14]), spx.norm.norms_l2_cw(x));
    verifyEqual(testCase, sqrt([3 12 27]'), spx.norm.norms_l2_rw(x));
    verifyEqual(testCase, [3 3 3], spx.norm.norms_linf_cw(x));
    verifyEqual(testCase, [1 2 3]', spx.norm.norms_linf_rw(x));
    verifyEqual(testCase, 18, spx.norm.cr_l1_l1(x));
    verifyEqual(testCase, 18, spx.norm.rc_l1_l1(x));
    verifyEqual(testCase, sqrt(108), spx.norm.cr_l1_l2(x));
    verifyEqual(testCase, sqrt(126), spx.norm.rc_l1_l2(x), 'AbsTol', 1e-6);
    verifyEqual(testCase, 6, spx.norm.cr_l1_linf(x));
    verifyEqual(testCase, 9, spx.norm.rc_l1_linf(x));
    verifyEqual(testCase, 3*sqrt(14),spx.norm.cr_l2_l1(x));
    verifyEqual(testCase, sqrt(3) + sqrt(12) + sqrt(27), spx.norm.rc_l2_l1(x));
    verifyEqual(testCase, sqrt(42), spx.norm.cr_l2_l2(x), 'AbsTol', 1e-6);
    verifyEqual(testCase, sqrt(42), spx.norm.rc_l2_l2(x), 'AbsTol', 1e-6);
    verifyEqual(testCase, sqrt(14), spx.norm.cr_l2_linf(x));
    verifyEqual(testCase, sqrt(27), spx.norm.rc_l2_linf(x));
    verifyEqual(testCase, 9, spx.norm.cr_linf_l1(x));
    verifyEqual(testCase, 6, spx.norm.rc_linf_l1(x));
    verifyEqual(testCase, sqrt(27), spx.norm.cr_linf_l2(x), 'AbsTol', 1e-6);
    verifyEqual(testCase, sqrt(14), spx.norm.rc_linf_l2(x), 'AbsTol', 1e-6);
    verifyEqual(testCase, 3, spx.norm.cr_linf_linf(x));
    verifyEqual(testCase, 3, spx.norm.rc_linf_linf(x));
end

author: indigits

repo: sparse-plex


text:
%HIST Display feature histograms
%
%     HANDLE = HIST(A,P,NX)
%     HANDLE = HIST(+A)
%
% INPUT
%  A       Dataset
%  P       # of bins
%  NX      # of histograms displayed in a row
%
% OUTPUT
%  HANDLE  handle of subplot
%
% DESCRIPTION
% For all feature (columns) of A a histogram is plot using P bins.
% These histograms are plot as subplots in a single figure, displaying
% NX histograms in a row. In HANDLE the handles of the subplots are
% returned.
%
% Note that this routine is not a true overload of the HIST command.
% Use HIST(+A) if that is desired.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% DEMSWISSROLLFULLLLE1 Demonstrate LLE on the oil data.

% MLTOOLS

[Y, lbls] = lvmLoadData('swissRollFull');

options = lleOptions(4, 2);
model = lleCreate(2, size(Y, 2), Y, options);
model = lleOptimise(model);

lvmScatterPlotColor(model, model.Y(:, 2));

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, model.Y(:, 2), 'SwissRollFull', 1, true);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [enhancedMag, enhancedWav] = RunLogMMSE(noisyFile, fs, isFile)
addpath('../../../../Enhancement/Loizou/MATLAB_code/statistical_based');

w = warning ('off','all');
if isFile
    noisyFile = ExtractWordsFromString_v2(noisyFile);
    noisyFile = noisyFile{1};
    % [enhancedWav, fs] = logmmse_vectorize(noisyFile);
    [enhancedWav, fs] = logmmse(noisyFile, [], isFile);
else
    [enhancedWav] = logmmse(noisyFile, fs, isFile);
end

if fs==16000
    fftlen = 512;
else
    fftlen = 256;
end

[~,enhancedMag] = wav2abs_multi(enhancedWav, fs, 0.01, 0.025, fftlen);
enhancedMag = 2*log(abs(enhancedMag(1:fftlen/2+1,:)));

end

author: singaxiong

repo: SignalGraph


text:
function besi0_test ( )

%*****************************************************************************80
%
%% BESI0_TEST tests R4_BESI0 and R8_BESI0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESI0_TEST:\n' );
  fprintf ( 1, '  Test BESI0_VALUES, R4_BESI0, R8_BESI0\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X       BESI0(X)\n' );
  fprintf ( 1, '                  R4_BESI0(X)         Diff\n' );
  fprintf ( 1, '                  R8_BESI0(X)         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = bessel_i0_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_besi0 ( single ( x ) );
    fx3 = r8_besi0 ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer


text:
classdef SettingsLevelSetCreatorFactory < handle

    methods (Access = public, Static)

        function obj = create(s)
            switch s.type
                case 'circle'
                    obj = SettingsLevelSetSphereNdim(s);
                case 'circleInclusion'
                    obj = SettingsLevelSetSphereNdim(s);
                case 'sphere'
                    obj = SettingsLevelSetSphereNdim(s);
                case 'sphereInclusion'
                    obj = SettingsLevelSetSphereNdim(s);
                case 'cylinder'
                    obj = SettingsLevelSetSphereNdim(s);
                case 'horizontal'
                    obj = SettingsLevelSetHorizontalInclusion(s);
                case 'squareInclusion'
                    obj = SettingsLevelSetSquareInclusion(s);
                case 'smoothSquare'
                    obj = SettingsLevelSetSmoothSquareInclusion(s);
                case {'rectangle', 'rectangleInclusion'}
                    obj = SettingsLevelSetRectangleInclusion(s);
                case 'smoothRectangle'
                    obj = SettingsLevelSetSmoothRectangleInclusion(s);
                case 'feasible'
                    obj = SettingsLevelSetCreator(s);
                case 'rand'
                    obj = SettingsLevelSetCreator(s);
                case 'holes'
                    obj = SettingsLevelSetSeveralHoles(s);
                case 'full'
                    obj = SettingsLevelSetCreator(s);
                case 'horizontalFibers'
                    obj = SettingsLevelSetHorizontalFibers(s);
                case 'given'
                    obj = SettingsLevelSetGiven(s);
                case 'Vigdergauz'
                    obj = SettingsLevelSetVigdergauz(s);
                otherwise
                    error('Invalid initial value of design variable.');
            end

        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan