text:
function [d] = triangle2distance(tri, pos, s, maxinitdist)

% TRIANGLE2DISTANCE computes the geodesic distance (across the edges) on a
% mesh, using Dijkstra's algorithm. The Dijkstra code is an efficient
% vectorized version of a function from MIT's graphtool toolbox, operating
% on an adjacency matrix.
%
% Use as
%   d = triangle2distance(tri, pos, s)
%
% Input arguments:
%   tri = Mx3 matrix describing the triangles
%   pos = Nx3 matrix describing the position of the vertices
%   s   = (can be empty), scalar or vector with indices for the points for
%         which the distance (to all other points) will be computed. If
%         empty or not defined, all points will be considered.
%
% Output argument:
%   d   = Nxnumel(s) distance matrix

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

adj = triangle2connectivity(tri, pos);
n   = length(adj);
if nargin<3 || isempty(s)
  s  = 1:n;
end
if nargin<4
  maxinitdist = inf;
end
ns = length(s);

d    = inf*ones(n,ns); % distance s-all nodes
for k = 1:ns
  d(s(k),k) = 0; % s-s distance
end

for k = 1:ns
  if mod(k,10)==0,
    fprintf('computing distance between node %d and nodes 1:%d\n', s(k), n);
  end
  T = 1:n;    % node set with shortest paths not found

  if isfinite(maxinitdist)
    initd = sqrt(sum((pos-pos(s(k),:)).^2,2));
    T(initd>=maxinitdist) = [];
  end

  %while ~isempty(T) %%% a for-loop goes ~10% faster, and we know how often
  %we need to iterate
  for m = 1:n
    if mod(m,1000)==0,
      fprintf('looping across vertices %d/%d\n', m, n);
    end
    [dmin,ind] = min(d(T,k));

    adj_ = adj(T,T(ind));
    d_   = d(T,k);

    % logic: shrink the distance if there's an edge (adj_>0) AND if it's
    % shorter to travel through this edge
    criterion = adj_ > 0 & d_ > d(T(ind),k)+adj_;
    d(T(criterion),k) = d(T(ind),k) + adj_(criterion);

    T(ind) = [];
  end
end

% the below code is from MIT's graphtool toolbox. the above code that
% computes the distance based on an adjacency matrix is taken from there,
% but the vectorized functionality is ~10 times as fast.

% Implements a simple version of the Dijkstra shortest path algorithm
% Returns the distance from a single vertex to all others, doesn't save the path
% INPUTS: adjacency matrix (adj), start node (s)
% OUTPUTS: shortest path length from start node to all other nodes
% Note: works with a weighted/directed matrix
% GB, Last Updated: December 13, 2004

% Copyright (c) 2011, Massachusetts Institute of Technology.
% All rights reserved.
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%
%    Redistributions of source code must retain the above copyright notice, this list
%    of conditions and the following disclaimer.
%    Redistributions in binary form must reproduce the above copyright notice, this list
%    of conditions and the following disclaimer in the documentation and/or other materials
%    provided with the distribution.
%    Neither the name of the Massachusetts Institute of Technology nor the names of its
%    contributors may be used to endorse or promote products derived from this software without
%    specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
% IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
% IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
% THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

% function d = simple_dijkstra(adj,s)
%
% n=length(adj);
% d = inf*ones(1,n); % distance s-all nodes
% d(s) = 0;    % s-s distance
% T = 1:n;    % node set with shortest paths not found
%
% while not(isempty(T))
%     [dmin,ind] = min(d(T));
%     for j=1:length(T)
%         if adj(T(ind),T(j))>0 & d(T(j))>d(T(ind))+adj(T(ind),T(j))
%             d(T(j))=d(T(ind))+adj(T(ind),T(j));
%         end
%     end
%     T = setdiff(T,T(ind));
%
% end

author: fieldtrip

repo: fieldtrip


text:
function [y,deriv] = shift_function(w,shift,f)
% This is an MV2DF (see MV2DF_API_DEFINITION.readme) which
% represents the new function,
%
%    g(w) = shift(w)+f(w),
%
% where shift is scalar-valued and f is matrix-valued.
%
%
% Here shift and f are function handles to MV2DF's.

if nargin==0
    test_this();
    return;
end

if isempty(w)

    s = stack(w,shift,f);
    map = @(s) s(2:end)+s(1);
    transmap = @(y) [sum(y);y];
    y = linTrans(s,map,transmap);
    return;
end

if isa(w,'function_handle')
    f = shift_function([],shift,f);
    y = compose_mv(f,w,[]);
    return;
end

f = shift_function([],shift,f);
if nargout==1
    y = f(w);
else
    [y,deriv] = f(w);
end

function test_this()

m = 5;
n = 10;
data = randn(m,n);
shift = 3;
w = [data(:);shift];

g = subvec([],m*n+1,1,m*n);
shift = subvec([],m*n+1,m*n+1,1);

f = shift_function([],shift,g);
test_MV2DF(f,w);

author: nesl


text:
function mne_license

% MNE_LICENSE prints the license only once upon the first call to
% this function. If the user does a "clear all", the license will
% again be shown.  This function should be included in every openmeeg
% function to ensure that the license is displayed at least once.

% Copyright (C) 2009-2011, Robert Oostenveld

persistent status

if isempty(status)
  % this is for the first time
  status = false;
end

if ~status
  % show the license
  try
    tmp = which('mne_license');
  catch
    tmp = mfilename;
  end
  [p, f] = fileparts(tmp);
  clc
  fprintf('==============================================================================\n');
  type(fullfile(p, 'LICENSE.txt'));
  fprintf('==============================================================================\n');
  % remember that it has been shown
  status = true;
  pause(3)
end

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function persist ( )

%*****************************************************************************80
%
%  PERSIST demonstrates the use of a peristent variable.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 January 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PERSIST\n' );
  fprintf ( 1, '  Demonstrate how a PERSISTENT variable can be used\n' );
  fprintf ( 1, '  in a lower level function, to store a value between calls.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I  Plus Total\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 5
    j = floor ( 10 * rand ( ) );
    total = sumup ( j );
    fprintf ( 1, '  %4d  %4d  %4d\n', i, j, total );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PERSIST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function total = sumup ( arg )

%*****************************************************************************80
%
%  SUMUP adds the input argument to an internal total and returns the new value.
%
%  Discussion:
%
%    The point of this function is simply to illustrate how a persistent
%    variable can be used to retain a value within a function between calls.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ARG, the value to be added to the total.
%
%    Output, real VALUE, the running total.
%
%  Local parameters:
%
%    Local, real VALUE, an internal persistent variable which records the
%    current running total.
%
  persistent value;

  if ( isempty ( value ) )
    value = 0;
  end

  value = value + arg;

  total = value;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% SVARGPLVMTOOLBOXES Load in the relevant toolboxes for shared variational gplvm.

importLatest('netlab');
importLatest('prior');
importLatest('optimi');
importLatest('datasets');
importLatest('mltools');
importLatest('kern');
importLatest('ndlutil');
importLatest('mocap');
importLatest('gp');
importLatest('noise');
importLatest('vargplvm');
importLatest('matlab');

author: SheffieldML

repo: vargplvm


text:
function sar = sar(high, low)
    %constants
    uptrend = 1;
    downtrend = 2;
    sar(1)=NaN;

    %variables
    alpha=0.02;
    EP=high(1);%between the high and low

    %'previous trade' is either long or short on day 1. sar(1)=NaN

    %decide whether to start with long or short
    if(high(1)>high(2))
        position=uptrend;
    else
        position=downtrend;
    end

    if(position==uptrend)
        sar(1)=high(1);
    else
        sar(1)=low(1);
    end

    for n=1:length(high)-1
        if(position==uptrend)
            EP=min(EP,low(n));
        else
            EP=max(EP,high(n));
        end
        alpha=min(alpha+0.02,0.2);
        sar(n+1)=sar(n)+alpha*(EP-sar(n));

        %"Parabolic SAR is never moved within the range of the current
        %or previous day (highest High to lowest Low over the 2 days)."
        if(n>2)
            if(sar(n+1)<max(high(n-1),high(n)) && sar(n)>min(low(n-1),low(n)))
                if(position==downtrend)
                    sar(n+1)=min(low(n-1),low(n));
                else
                    sar(n+1)=max(high(n-1),high(n));
                end
            end
        end
        if(sar(n+1)<high(n+1) && sar(n+1)>low(n+1))
            sar(n+1)=EP;
            if(position==uptrend)
                position=downtrend;
                EP=low(n);
            else
                position=uptrend;
                EP=high(n);
            end
            alpha=0.02;
            %sar(n+1)=sar(n)+alpha*(EP-sar(n));
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function [tParam, regParam] = histology_intraframe_reg(pathstr, files, opts)
% histology_intraframe_reg  Registration diffusion of histology slices.
%
% histology_intraframe_reg computes the registration of each slices to the
% two adjacent slices. The two resulting transforms are combined to
% produce a compromise transform. The process is applied iteratively.
%
% [TPARAM, REGPARAM] = histology_intraframe_reg(PATHSTR, FILES)
%
%   PATHSTR is the full path to the files. If empty, PATHSTR='.'.
%
%   FILES is the result of a dir() command, e.g. dir('*.png'). The
%   function expects a list of histology images.
%
%   TPARAM is an array of structs with the transforms for each image, where
%   TPARAM(I) is the nested transform for the I-th image. Note that if
%   ta is a series of nested transforms, e.g.
%
%     ta.Transform = 'EulerTransform';
%     ta.InitialTransformParametersFileName = tb;
%
%     tb.Transform = 'BSplineTransform';
%     tb.InitialTransformParametersFileName = 'NoInitialTransform';
%
%   tb is the "initial transform" of ta, but because ITK and elastix define
%   the transform of an image in the inverse direction, in practice what
%   happens is that ta is applied *first* to the image, and then tb is
%   applied to the result.
%
% ... = histology_intraframe_reg(..., OPTS)
%
%   OPTS is a struct with algorithm parameters:
%
%     'MaxIter':   (def 10) Maximum number of diffusion iterations.
%
%     'DiffusionCoefficient': (def 0.4) Value of diffusion coefficient.
%
%     'verbose':   (default 0) Show elastix output on the screen.
%
%     't0':        (default []) Initial transform in elastix format for the
%                  images (typically rigid transforms to provide an initial
%                  rough alignment).
%
% See also: histology_preprocessing.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2014 University of Oxford
% Version: 0.2.2
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 3);
nargoutchk(0, 2);

if (length(files) < 3)
    error('At least three files are needed for registration')
end

% defaults
if (nargin < 3 || isempty(opts) || ~isfield(opts, 'MaxIter'))
    opts.MaxIter = 10;
end
if (isempty(opts) || ~isfield(opts, 'DiffusionCoefficient'))
    opts.DiffusionCoefficient = 0.4;
end
if (isempty(opts) || ~isfield(opts, 'verbose'))
    opts.verbose = 0;
end
if (isempty(opts) || ~isfield(opts, 't0'))
    opts.t0 = [];
end
if (isempty(pathstr))
    pathstr = '.';
end

if (~isempty(opts.t0) && (length(opts.t0) ~= length(files)))
    error('If opts.t0 initial transforms are provided, there must be one per image')
end

% create a struct with the registration parameters
regParam = generate_registration_parameters();

% create a duplicate of the images that we will deform repeteadly in the
% diffusion process
pathtemp = tempname;
mkdir(pathtemp);
for I = 1:length(files)
    copyfile([pathstr filesep files(I).name], ...
        [pathtemp filesep files(I).name]);
end

% diffusion computation loop
for S = 1:opts.MaxIter

    % initial transform to apply to images before registration in this
    % iteration
    if (S==1)

        t0 = opts.t0;

    else

        t0 = tParam;

    end

    % compute registration diffusion for the first slice of the histology
    % volume
    tic
    disp(['S = ' num2str(S) '/' num2str(opts.MaxIter) ...
        ', I = ' num2str(1) '/' num2str(length(files))])
    tParam(1) = process_extreme_slices(...
        [pathstr filesep files(2).name], ... % fixed
        [pathstr filesep files(1).name], ... % moving
        t0(2), t0(1), regParam, opts);
    disp(['time = ' num2str(toc) 'sec'])

    % loop for all the other slices
    for I = 2:length(files)-1

        tic
        disp(['S = ' num2str(S) '/' num2str(opts.MaxIter) ...
            ', I = ' num2str(I) '/' num2str(length(files))])

        tParam(I) = process_intermediate_slices(...
            [pathstr filesep files(I-1).name], ... % fixed a
            [pathstr filesep files(I).name], ...   % moving
            [pathstr filesep files(I+1).name], ... % fixed b
            t0(I-1), t0(I), t0(I+1), regParam, opts);
        disp(['time = ' num2str(toc) 'sec'])

    end

    % compute registration diffusion for the last slice of the histology
    % volume
    N = length(files);
    tic
    disp(['S = ' num2str(S) '/' num2str(opts.MaxIter) ...
        ', I = ' num2str(length(files)) '/' num2str(length(files))])
    tParam(N) = process_extreme_slices(...
        [pathstr filesep files(N-1).name], ... % fixed
        [pathstr filesep files(N).name], ...   % moving
        t0(N-1), t0(N), regParam, opts);
    disp(['time = ' num2str(toc) 'sec'])

end

% delete temp files
rmdir(pathtemp, 's')

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate_registration_parameters()
%
%   Create a struct with the parameters necessary for registration
function regParam = generate_registration_parameters()

regParam.FixedInternalImagePixelType = 'float';
regParam.FixedImageDimension = 2;
regParam.MovingInternalImagePixelType = 'float';
regParam.MovingImageDimension = 2;
regParam.Registration = 'MultiMetricMultiResolutionRegistration';
regParam.FixedImagePyramid = {'FixedSmoothingImagePyramid', ...
    'FixedSmoothingImagePyramid', 'FixedSmoothingImagePyramid'};
regParam.MovingImagePyramid = {'MovingSmoothingImagePyramid', ...
    'MovingSmoothingImagePyramid', 'MovingSmoothingImagePyramid'};
regParam.Interpolator = {'BSplineInterpolator', 'BSplineInterpolator', ...
    'BSplineInterpolator'};
regParam.Metric = {'AdvancedMeanSquares', 'AdvancedMeanSquares', ...
    'AdvancedMeanSquares'};
regParam.Optimizer = 'RegularStepGradientDescent';
regParam.ResampleInterpolator = 'FinalBSplineInterpolator';
regParam.Resampler = 'DefaultResampler';
regParam.Transform = 'BSplineTransform';
regParam.BSplineTransformSplineOrder = 3;
regParam.UseCyclicTransform = 'false';
regParam.AutomaticTransformInitialization = 'false';
regParam.AutomaticTransformInitializationMethod = 'GeometricalCenter';
regParam.UseFastAndLowMemoryVersion = 'true';
regParam.UseJacobianPreconditioning = 'false';
regParam.FiniteDifferenceDerivative = 'false';
regParam.MaximumNumberOfSamplingAttempts = 10;
regParam.ErodeMask = 'false';
regParam.NumberOfResolutions = 3;
regParam.ImagePyramidSchedule = [16 16 4 4 1 1];
regParam.FinalGridSpacingInVoxels = 32;
regParam.GridSpacingSchedule = [16 4 1];
regParam.HowToCombineTransforms = 'Compose';
regParam.UseDirectionCosines = 'true';
regParam.MaximumNumberOfIterations = 400;
regParam.MaximumStepLength = 50;
regParam.MinimumStepLength = 0.500;
regParam.MinimumGradientMagnitude = 1e-3;
regParam.NumberOfFixedHistogramBins = 32;
regParam.NumberOfMovingHistogramBins = 32;
regParam.NumberOfHistogramBins = 32;
regParam.FixedLimitRangeRatio = 0;
regParam.MovingLimitRangeRatio = 0;
regParam.FixedKernelBSplineOrder = 3;
regParam.MovingKernelBSplineOrder = 3;
% regParam.ImageSampler = {'RandomSparseMask' 'RandomSparseMask' 'RandomSparseMask'};
regParam.ImageSampler = {'RandomCoordinate' 'RandomCoordinate' 'RandomCoordinate'};
regParam.NewSamplesEveryIteration = 'true';
regParam.NumberOfSpatialSamples = [0 0 0]; % to be set for each slice
regParam.FixedImageBSplineInterpolationOrder = 1;
regParam.UseRandomSampleRegion = 'false';
regParam.BSplineInterpolationOrder = 1;
regParam.FinalBSplineInterpolationOrder = 1;
regParam.DefaultPixelValue = 0;
regParam.WriteResultImage = 'true';
regParam.WriteResultImageAfterEachResolution = 'false';
regParam.WriteTransformParametersEachIteration = 'false';
regParam.WriteTransformParametersEachResolution = 'true';
regParam.ResultImagePixelType = 'unsigned char';
regParam.ResultImageFormat = 'png';
regParam.CompressResultImage = 'true';

end

% process_extreme_slices()
%
%   The first and last slices of the histology volume don't have slices on
%   both sides. This function deals with those two cases (it has to be run
%   twice, once for the first slice, and once for the last slice).
%
%   *M refers to the moving slice, i.e. the first or last histology slice.
%   *F refers to the fixed adjacent slice.
%
%   FILEM, FILEF is the path and filename to the histology images.
%
%   T0M, T0F is the nested (accumulated) transform to apply to each image
%   before registering them.
%
%   REGPARAM is the struct with the registration parameters for elastix.
%
%   OPTS is a struct with optional algorithm parameters.
function t0m = process_extreme_slices(filef, filem, t0f, t0m, regParam, opts)

DEBUG = 0;

% load images
imm = imread(filem);
imf = imread(filef);

% preprocess histograms
[imm, imf] = histology_preprocessing(imm, imf);

% apply accumulated initial transform (make sure that added voxels are
% black, so that they blend in with the background)
t0m.DefaultPixelValue = 0;
t0f.DefaultPixelValue = 0;
imm = transformix(t0m, imm);
imf = transformix(t0f, imf);

% create masks
opts.mMask = uint8(sum(imm, 3) > 0);
opts.fMask = uint8(sum(imf, 3) > 0);

% remove some background noise in the mask
se = strel('disk', 1);
opts.mMask = imerode(opts.mMask, se);
opts.fMask = imerode(opts.fMask, se);

% dilate mask to cover some background
se = strel('disk', 10);
opts.mMask = imdilate(opts.mMask, se);
opts.fMask = imdilate(opts.fMask, se);

% we use 10% of the voxels within the mask
nSamples = nnz(imf) / 10 * regParam.GridSpacingSchedule(end);
regParam.NumberOfSpatialSamples = round(nSamples ./ regParam.GridSpacingSchedule);

% register slices
[tma, immreg, iterinfo] = elastix(regParam, imf, imm, opts);

if (DEBUG)
    subplot(2, 1, 1)
    hold off
    imshowpair(imf, imm)
    subplot(2, 1, 2)
    hold off
    imshowpair(imf, immreg)
end

% create a new combined transform as a linear combination of the two
% computed transforms, and apply diffusion coefficient (the new transform
% has to go at the back of the list so that it is applied last)
tma.TransformParameters ...
    = 2 * tma.TransformParameters * opts.DiffusionCoefficient;
str = 't0m';
while (isstruct(eval([str '.InitialTransformParametersFileName'])))
    str = [str '.InitialTransformParametersFileName'];
end
eval([str '.InitialTransformParametersFileName = tma;']);

end

% process_intermediate_slices()
%
%   All slices of the histology volume except the first and last have
%   slices on both sides. This function computes the difussion registration
%   of a slice to its two adjacent slices.
%
%   *M refers to the moving slice, i.e. the first or last histology slice.
%   *FA, *FB refer to the two fixed adjacent slices.
%
%   FILEM, FILEFA, FILEFB are the path and filename to the histology
%   images.
%
%   T0M, T0FA, T0FB are the nested (accumulated) transforms to apply to
%   each image before registering them.
%
%   REGPARAM is the struct with the registration parameters for elastix.
%
%   OPTS is a struct with optional algorithm parameters.
function t0m = process_intermediate_slices(filefa, filem, filefb, t0fa, t0m, t0fb, regParam, opts)

DEBUG = 0;

% load images
imm = imread(filem);
imfa = imread(filefa);
imfb = imread(filefb);

% preprocess histograms
[~, imfa] = histology_preprocessing(imm, imfa);
[imm, imfb] = histology_preprocessing(imm, imfb);

% apply accumulated initial transform (make sure that added voxels are
% black, so that they blend in with the background)
t0m.DefaultPixelValue = 0;
t0fa.DefaultPixelValue = 0;
t0fb.DefaultPixelValue = 0;
imm = transformix(t0m, imm);
imfa = transformix(t0fa, imfa);
imfb = transformix(t0fb, imfb);

% create masks
opts.mMask = create_mask(imm);
fMaska = create_mask(imfa);
fMaskb = create_mask(imfb);

% register slices
opts.fMask = fMaska;
regParam.NumberOfSpatialSamples = round(nnz(fMaska) / 10); % use 10% of the voxels within the mask
[tma, immrega, iterinfoa] = elastix(regParam, imfa, imm, opts);
opts.fMask = fMaskb;
regParam.NumberOfSpatialSamples = round(nnz(fMaskb) / 10); % use 10% of the voxels within the mask
[tmb, immregb, iterinfob] = elastix(regParam, imfb, imm, opts);

% visualize registered images
if (DEBUG)
    subplot(2, 2, 1)
    hold off
    imshowpair(imfa, imm)
    subplot(2, 2, 3)
    hold off
    imshowpair(imfa, immrega)
    subplot(2, 2, 2)
    hold off
    imshowpair(imfb, imm)
    subplot(2, 2, 4)
    hold off
    imshowpair(imfb, immregb)
end

% create a new combined transform as a linear combination of the two
% computed transforms, and apply diffusion coefficient (the new transform
% has to go at the back of the list so that it is applied last)
tma.TransformParameters ...
    = (tma.TransformParameters + tmb.TransformParameters) ...
    * opts.DiffusionCoefficient;
str = 't0m';
while (isstruct(eval([str '.InitialTransformParametersFileName'])))
    str = [str '.InitialTransformParametersFileName'];
end
eval([str '.InitialTransformParametersFileName = tma;']);

% visualize deformation field
if (DEBUG)

    % create grid image
    aux = zeros(size(imm), 'uint8');
    for I = 1:5
        aux(I:100:end, :, :) = 255;
        aux(:, I:100:end, :) = 255;
    end

    % deform grid according to both transforms independently, and the
    % combined one
    auxrega = transformix(tma, aux);
    auxregb = transformix(tmb, aux);
    auxregm = transformix(tm, aux);

    % display 3 deformation fields
    subplot(3, 1, 1)
    imagesc(auxrega)
    axis equal
    subplot(3, 1, 2)
    imagesc(auxregb)
    axis equal
    subplot(3, 1, 3)
    imagesc(auxregm)
    axis equal

end

end

function mask = create_mask(im)

% threshold image
mask = uint8(sum(im, 3) > 0);

% detect connected components of the mask
cc = bwconncomp(mask);

% remove small components
stats = regionprops(cc, 'Area');
idx = [stats.Area] < 5000;
mask(cat(1, cc.PixelIdxList{idx})) = 0;

% smooth mask and dilate a bit
se = strel('disk', 15);
mask = imdilate(mask, se);
se = strel('disk', 5);
mask = imerode(mask, se);

% fill in the holes
mask = imfill(mask, 'holes');

% make sure we only have one connected components
mask = bwrmsmallcomp(mask);

% fit ellipse to the segmentation mask
stats = regionprops(mask, 'MajorAxisLength', 'MinorAxisLength', ...
    'Centroid', 'Orientation');

% simplify notation
xc = stats.Centroid(1);
yc = stats.Centroid(2);
a = stats.MajorAxisLength / 2;
b = stats.MinorAxisLength / 2;
alpha = stats.Orientation;

% create elliptical segmentation mask
[x, y] = meshgrid(1:size(mask, 2), 1:size(mask, 1));
xy = [x(:) - xc, y(:) - yc] ...
    * [cosd(alpha) sind(alpha); -sind(alpha) cosd(alpha)];

mask = (b^2 * xy(:, 1).^2 + a^2 * xy(:, 2).^2 <= (a * b)^2);
mask = reshape(mask, size(im, 1), size(im, 2));

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
  function [xs, wlsden] = wls_gcd(A, W, yy, xhat, niter, ng, chat)
%|function [xs, wlsden] = wls_gcd(A, W, yy, xhat, niter, ng, chat)
%| weighted least squares minimized by
%| grouped coordinate descent algorithm
%| in
%|	A	[nn np]			system matrix
%|	W	[nn nn]			inverse data covariance
%|	yy	[nn 1]			noisy data
%|	nder1	[nn 1]
%|	xhat	[np 1] or [nx ny]	initial estimate
%|	niter		# total iterations
%|	ng		# groups, or, if array, groups
%| out
%|	xs	[np niter]		estimates each iteration
%|
%| Copyright July 1996, Jeff Fessler, University of Michigan

if nargin < 6, help(mfilename), error(' '), end

if 1 ~= exist('chat'), chat = 1; end
%chat = (1==exist('mask'));

[nx,ny] = size(xhat);

if 1~=exist('mask') || isequal(mask, 1)
	mask = ones(nx,ny);
end

np = ncol(A);

if numel(xhat) ~= np
	xhat = xhat(mask(:));
end

if ~isempty(nder1)
	back_nder1 = A' * nder1(:);
else
	back_nder1 = zeros(np,1);
end
yy = yy(:);
xs = zeros(np, niter);
xs(:,1) = xhat;

if chat
	im(321, embed(xhat,mask), 'x0')
end

% build groups
if numel(ng) == 1
	groups	= zeros(np, ng);
	for ig = 1:ng
		gg = [ig:ng:np]';
		groups(gg,ig) = ones(size(gg));
	end
elseif numel(ng) == 2
	groups	= group2d(nx, ny, ng(1), ng(2));
	ng	= ncol(groups);
else
	groups	= ng;
	ng	= ncol(groups);
end

% overlapping groups?
if any(sum(groups')) > 1, error 'overlap groups', end

if chat
	im(322, embed(groups * [1:ng]',mask), 'Groups')
end

% precompute quadratic part of denominator
% change '1' to 'ig' for overlapping groups
wlsden	= zeros(np, 1);
if 1
	for ig = 1:ng
		gg = groups(:,ig);
		printm('precompute %g, assuming nonnegative A', ig)

%		wlsden(gg,1) = abs(A(:,gg))' * (W * sum(abs(A(:,gg)'),1)');
%		wlsden(gg,1) = A(:,gg)' * (W * sum(A(:,gg)')',1);

		t = A * gg;	% faster than: sum(A(:,gg)',1)';
		t = W * t;	% only slightly slower than: diag(W) .* t;
		%	faster than either: A(:,gg)'*t or A' * t !
		t = t' * A;
		wlsden(gg,1) = t(:,gg)';
	end
else
	wlsden	= ones(np, 1);	% for testing
end

if chat
	im(323, embed(wlsden,mask), 'WLS Denominator')
end

% compare GCA and SCA denominators (for WLS part)
if chat && (np < 1000)
	t = full(sum(A .* (W * A))'); % SCA den
	t = t(wlsden ~= 0) ./ wlsden(wlsden ~= 0);
	if min(t(:)) < 1-10*eps
		printm(['Under-relax. range:' sprintf(' %g', minmax(t))])
		t = embed(t,wlsden);
		im(324, embed(t,mask), 'Under-relaxer')
	else
		printm 'Note: GCA with ideal denominator!'
	end
end

% 'depierro' factors for penalty? (GCA vs SCA)
for ig = 1:ng
	gg = groups(:,ig);
	t = abs(C(:,gg))' * sum(abs(C(:,gg)'),1)';
	t = sum(C(:,gg).^2)' ./ t; % SCA / GCA
	t = minmax(t);
	if min(t) < 1-10*eps
		printm(['Virtual under-relax. range:' sprintf(' %g', t)])
		error 'depierro not done'
	end
	clear t
end

if chat
	im(325, embed(sum(C.^2)',mask), 'sum(C.^2)')
	input('hit key');
end

res = A * xhat - yy;

for ii = 2:niter

	for ig = 1:ng
		gg = groups(:,ig);

		xn = xhat(gg);
		xw = xn;

		ader = ( (W * res)' * A(:,gg) )' + back_nder1(gg);

		for jj = 1:nsub
			dx = C * xhat;
			eval(sprintf(wstring, 'wt', 'dx'));
			denom = wlsden(gg) + (C(:,gg).^2)' * wt;
			num = ader + wlsden(gg).*(xw-xn) ...
				+ C(:,gg)' * (wt .* dx);
			xw = xw - num ./ denom;
			xhat(gg) = xw;
		end

		res = res + A(:,gg) * (xw - xn);
	end

	xs(:,ii) = xhat;
	if chat && (ii < 10 || rem(ii,10) == 0)
		im(326, embed(xhat,mask), 'x hat')
		if ny == 1
			im(324, wt, 'wts')
		end
	end
	if chat
		printm('max %g %g', max(xhat), max(abs(xhat-xs(:,ii-1)))/max(xhat)*100)
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function cdf = genlogistic_cdf ( x, a, b, c )

%*****************************************************************************80
%
%% GENLOGISTIC_CDF evaluates the Generalized Logistic CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < B,
%    0.0 < C.
%
%    Output, real CDF, the value of the CDF.
%
  y = ( x - a ) / b;

  cdf = 1.0 / ( 1.0 + exp ( - y ) )^c;

  return
end


text:
function quadrule_test087 ( )

%*****************************************************************************80
%
%% QUADRULE_TEST087 tests HERMITE_EK_COMPUTE.
%
%  Discussion:
%
%    I used this test to generate tabular values of weights and
%    abscissas for Gauss-Hermite quadrature.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2011
%
%  Author:
%
%    John Burkardt
%
  n = 31;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADRULE_TEST087\n' );
  fprintf ( 1, '  HERMITE_EK_COMPUTE computes Gauss-Hermite data;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute the data for N = %d\n', n );

  [ x, w ] = hermite_ek_compute ( n );

  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '    x(%d) = %24.16e;\n', i, x(i) );
  end
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '    w(%d) = %24.16e;\n', i, w(i) );
  end

  return
end


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [glob_num, glob_ndof] = mp_interface_3d (interfaces, sp)

warning ('geopdes:obsolete', 'Function MP_INTERFACE_3D is obsolete. Using MP_INTERFACE instead')

[glob_num, glob_ndof] = mp_interface (interfaces, sp);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function []=strsdisp(decimaldates1,stringdates1,strshocks_estimates,endo,pref,IRFt,strctident)

if IRFt==1||IRFt==2||IRFt==3
    identified=size(endo,1); % fully identified
    labels=endo; % simply use the name of the endogenous variables as labels
elseif IRFt==4 || IRFt==6 %if the model is identified by sign restrictions or sign restrictions (+ IV)
    identified=size(strctident.signreslabels_shocks,1); % count the labels provided in the sign res sheet (+ IV)
    labels=strctident.signreslabels_shocks; % signreslabels, only identified shocks
elseif IRFt==5
    identified=1; % one IV shock
    labels{identified,1}=strcat('IV Shock (',strctident.Instrument,')'); % one label for the IV shock
end

if pref.plot==1
% create structural shocks figure
strshocks=figure('Tag','BEARresults');
set(strshocks,'Color',[0.9 0.9 0.9]);
set(strshocks,'name','structural shocks');
ncolumns=ceil(identified^0.5);
nrows=ceil(identified/ncolumns);
for ii=1:identified
subplot(nrows,ncolumns,ii);
hold on
Xpatch=[decimaldates1' fliplr(decimaldates1')];
Ypatch=[strshocks_estimates{ii,1}(1,:) fliplr(strshocks_estimates{ii,1}(3,:))];
Fpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
set(Fpatch,'facealpha',0.6);
set(Fpatch,'edgecolor','none');
plot(decimaldates1,strshocks_estimates{ii,1}(2,:),'Color',[0.4 0.4 1],'LineWidth',2);
plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
hold off
minband=min(strshocks_estimates{ii,1}(1,:));
maxband=max(strshocks_estimates{ii,1}(3,:));
space=maxband-minband;
if space==0
    space=0.001;
end
Ymin=minband-0.2*space;
Ymax=maxband+0.2*space;
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'YLim',[Ymin,Ymax],'FontName','Times New Roman');
% create title
title(labels{ii,1},'FontName','Times New Roman','FontSize',10,'FontWeight','normal','interpreter','none');
end

end % pref.plot

% finally, record the results on excel
if pref.results==1
bear.data.excelrecord9fcn(stringdates1, labels, strshocks_estimates, identified, pref)
end

author: european-central-bank

repo: BEAR-toolbox


text:
function calpak_test69 ( )

%*****************************************************************************80
%
%% CALPAK_TEST69 tests YMDF_DIF_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2013
%
%  Author:
%
%    John Burkardt
%
  dhi = 1;
  dlo = 1;
  fhi = 0.0;
  flo = 0.0;
  mhi = 1;
  mlo = 1;
  seed = 123456789;
  yhi = 1960;
  ylo = 1970;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST69\n' );
  fprintf ( 1, '  YMDF_DIF_COMMON gets the day difference \n' );
  fprintf ( 1, '  between YMDF dates.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  YMDF1        YMDF2        (YMDF2 - YMDF1)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ y1, m1, d1, f1, seed ] = ymdf_uniform_common ( ylo, mlo, dlo, flo, ...
      yhi, mhi, dhi, fhi, seed );

    s1 = ymdf_to_s_common ( y1, m1, d1, f1 );

    [ y2, m2, d2, f2, seed ] = ymdf_uniform_common ( ylo, mlo, dlo, flo, ...
      yhi, mhi, dhi, fhi, seed );

    s2 = ymdf_to_s_common ( y2, m2, d2, f2 );

    [ days, ierror ] = ymdf_dif_common ( y1, m1, d1, f1, y2, m2, d2, f2 );

    fprintf ( 1, '  %10s  %10s  %11.2f\n', s1, s2, days );

  end

  return
end


text:
function value = hankel_n_condition ( n )

%*****************************************************************************80
%
%% HANKEL_N_CONDITION returns the L1 condition of the HANKEL_N matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real VALUE, the L1 condition.
%
  v = zeros ( n, 1 );

  v(1) = 1.0 / n;
  for i = 2 : n
    v(i) = 0.0;
    for j = 1 : i - 1
      v(i) = v(i) - ( n + j - i ) * v(j);
    end
    v(i) = v(i) / n;
  end

  a_norm = ( n * ( n + 1 ) ) / 2;
  b_norm = sum ( abs ( v(1:n) ) );
  value = a_norm * b_norm;

  return
end


text:
function f = gaussfunction(x,p)
% f = gaussfunction(x,p)
% p(1) = mean
% p(2) = std deviation (sqrt(var))
% p(3) = scale coefficient (if p(3)==1 then it is normalised)
normconst = 1/(sqrt(2*pi*p(2)^2));
f = p(3)*normconst*exp(-((x-p(1)).^2)/(2*p(2)^2));

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function model = yalmip2nonlinearsolver(model)

global newmodel

newmodel = 1;

model.dense = 0;
if ~(model.equalitypresolved || model.presolved)
    model = propagate_bounds_from_equalities(model);
end

K = model.K;
lb = model.lb;
ub = model.ub;
x0 = model.x0;
c = model.c;

% Pick out the positive conditions from cones ||Ax+b|| <= c'*x+d which will
% be treated as (Ax+b)'*(ax+b) <= (c'*x+d)^2,  c'*x+d >= 0
model = bounds_from_cones_to_lp(model);

if isempty(model.evaluation_scheme)
    model = build_recursive_scheme(model);
end
model = compress_evaluation_scheme(model);

% Do some pre-calc to be used in calls from fmincon
nonlinearindicies = union(find(model.variabletype~=0),model.evalVariables);
linearindicies    = setdiff(find(model.variabletype==0),nonlinearindicies);
model.nonlinearindicies = nonlinearindicies;
model.linearindicies    = linearindicies;

model.Anonlinineq = [];
model.bnonlinineq = [];
model.Anonlineq = [];
model.bnonlineq = [];

% Extract linear and nonlinear equality constraints
if K.f>0
    Aeq = -model.F_struc(1:1:K.f,2:end);
    beq = model.F_struc(1:1:model.K.f,1);

    nonlinear_equalities_indicies = find(any(Aeq(:,nonlinearindicies),2));
    model.Anonlineq = Aeq(nonlinear_equalities_indicies,:);
    model.bnonlineq = beq(nonlinear_equalities_indicies);

    Aeq(nonlinear_equalities_indicies,:) = [];
    beq(nonlinear_equalities_indicies,:) = [];
    Aeq(:,nonlinearindicies) = [];
    model.F_struc(1:model.K.f,:) = [];
    model.K.f = 0;
else
    Aeq = [];
    beq = [];
end

% Find nonlinear eualities implied by lower and upper bounds
if ~isempty(ub) && ~isempty(lb)
    nonlinearequality = find(lb(nonlinearindicies) == ub(nonlinearindicies));
    if ~isempty(nonlinearequality)
        for i = 1:length(nonlinearequality)
          %  model.Anonlineq = [model.Anonlineq;eyev(length(c),nonlinearindicies(nonlinearequality(i)))'];
            model.Anonlineq = [model.Anonlineq;sparse(1,nonlinearindicies(nonlinearequality(i)),1,1,length(c))];
            model.bnonlineq = [model.bnonlineq;lb(nonlinearindicies(nonlinearequality(i)))];
        end
    end
end

% Extract linear and nonlinear inequality constraints
if any(model.K.l)
    A = -model.F_struc(1:model.K.l,2:end);
    b = model.F_struc(1:model.K.l,1);

    nonlinear_inequalities_indicies = find(any(A(:,nonlinearindicies),2));

    model.Anonlinineq = A(nonlinear_inequalities_indicies,:);
    model.bnonlinineq = b(nonlinear_inequalities_indicies);

    A(nonlinear_inequalities_indicies,:) = [];
    b(nonlinear_inequalities_indicies,:) = [];
    A(:,nonlinearindicies) = [];

    model.F_struc(1:model.K.l,:) = [];
    model.K.l = 0;
else
    A = [];
    b = [];
end

% This helps with robustness in bnb in some cases
x0candidate = zeros(length(c),1);
if ~isempty(lb) && ~isempty(ub)
    bounded = find(~isinf(lb) & ~isinf(ub));
    x0candidate(bounded) = (lb(bounded) + ub(bounded))/2;
    bounded_below = find(~isinf(lb) & isinf(ub));
    x0candidate(bounded_below) = lb(bounded_below) + 0.5;
    bounded_above = find(~isinf(lb) & isinf(ub));
    x0candidate(bounded_above) = lb(bounded_above) + 0.5;
end

if isempty(x0)
    x0 = x0candidate(linearindicies);
else
    if ~isempty(lb) && ~isempty(ub)
        x0((x0 < lb) | (x0 > ub)) = x0candidate((x0 < lb) | (x0 > ub));
    end
    x0 = x0(linearindicies);
end

if ~isempty(lb)
    lb = lb(linearindicies);
end
if ~isempty(ub)
    ub = ub(linearindicies);
end

lb_old = lb;
ub_old = ub;
[lb,ub,A,b] = remove_bounds_from_Ab(A,b,lb,ub);
[lb,ub,Aeq,beq] = remove_bounds_from_Aeqbeq(Aeq,beq,lb,ub);

if any(model.variabletype == 4)
    problematic = find(any(model.monomtable(:,linearindicies) < 0 ,1));
    if ~isempty(problematic)
        problematic = problematic(find(x0(problematic)==0));
        Oneisfeas = problematic(find(ub(problematic) > 1));
        x0(Oneisfeas) = 1;
    end

    problematic = find(any(model.monomtable(:,linearindicies)~=fix(model.monomtable(:,linearindicies)) ,1));
    lb(problematic) = max(lb(problematic),0);
end
x0(find(lb==ub)) = lb(find(lb==ub));

if size(A,1) == 0
    A = [];
end

if size(b,1) == 0
    b = [];
end

if size(Aeq,1) == 0
    Aeq = [];
end

if size(beq,1) == 0
    beq = [];
end

if model.presolveequalities
    if ~isempty(beq) &  (~model.equalitypresolved | ~(isequal(lb,lb_old) & isequal(ub,ub_old)))
        % This helps when there are artificial variables introduced to model
        % nonlinear operators such as log(2*x+1)
        p.F_struc = [beq -Aeq];
        p.K.f = size(beq,1);
        p.lb = lb;
        p.ub = ub;
        p.variabletype = zeros(1,length(lb));
        p.binary_variables = [];
        p.integer_variables = [];
        p = propagate_bounds_from_equalities(p);
        lb = p.lb;
        ub = p.ub;
    end
end

model.A = A;
model.b = b;
model.Aeq = Aeq;
model.beq = beq;
model.lb = lb;
model.ub = ub;
model.x0 = x0;

model = setup_fmincon_params(model);

% Check if all derivatives are available
model.derivative_available = 1;
for i = 1:length(model.evalMap)
    if isempty(model.evalMap{i}.properties.derivative)
        model.derivative_available = 0;
        break
    end
end

% Some precomputation of computational scheme for Jacobian
allA = [model.Anonlineq;model.Anonlinineq];
if anyCones(model.K)
    allA = [allA;model.F_struc(startofSOCPCone(model.K):end,2:end)];
end
requested = any(allA',2);
[i,j] = find((model.deppattern(find(requested),:)));
requested(j) = 1;
if ~isempty(model.evalMap)
    % Recursive stuff is only possible if we have evaluation-based
    % operators
    model.Crecursivederivativeprecompute = precomputeDerivative(model,requested);
end

% Some precomputation of computational scheme for gradient
requested = model.c | any(model.Q,2);
[i,j,k] = find((model.deppattern(find(requested),:)));
requested(j) = 1;
model.frecursivederivativeprecompute = precomputeDerivative(model,requested);

% Precomputed list of bilinear expressions, used in
% apply_recursive_differentiation
model = compile_bilinearslist(model);
model = compile_quadraticslist(model);

model.binary_variables  = find(ismember(linearindicies,model.binary_variables));
model.integer_variables  = find(ismember(linearindicies,model.integer_variables));
model.semicont_variables  = find(ismember(linearindicies,model.semicont_variables));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% op_CSIoverlayMRI.m
% Brenden Kadota, Sunnybrook 2022
%
% DESCRIPTION:
% Overlay MRSI structure on spm_image graphic. Displays the MRSI specs onto
% the axial plane of the brain.
%
% USAGE:
% op_CSIoverlayMRI('/path/to/nifit/mri.nii', MRSIStruct, 1);
%
%
% INPUT:
% mriFileName           = string path to nifti file
% in                    = MRSI structure to overlay. In FID-A format
% coilIndex             = coilIndex to plot if coils have not been combined (default value = 1)
% averageIndex          = average index to plot
%
% OUTPUT:
%

function op_CSIoverlayMRI(MRSIStruct, mriFileName, coilIndex, averageIndex)
    arguments
        MRSIStruct (1, 1) struct
        mriFileName (:, 1) string
        coilIndex (1, 1) double = 1
        averageIndex (1, 1) double = 1
    end

    checkArguments(MRSIStruct);
    %call spm global variable
    global st
    global voxels
    global showVoxelOutlineFlag

    showVoxelOutlineFlag = false;

    % ensure filename is a character array not string
    mriFileName = char(mriFileName);
    if(contains(mriFileName, '~'))
        mriFileName = fullfile(getuserdir, mriFileName(2:end));
    end
    % display image using spm
    spm_image('display', mriFileName)

    % adds MRSI overlay function to SPM callback function
    st = addMRSIOverlayToSPMCallback(st, 'overlaySpectraOnMRI()');

    spmFigure = st.fig;
    buildMouseClickDropdown(spmFigure);
    buildPPMBoundTextBox(spmFigure);
    buildPlotTypeDropdown(spmFigure);
    buildShowVoxelButton(spmFigure);

    voxels = buildArrayOfVoxels(MRSIStruct);
    %plot the MRSI onto the spm mri graphic along the axial dimension.
    overlaySpectraOnMRI('real', min(MRSIStruct.ppm), max(MRSIStruct.ppm), MRSIStruct, coilIndex, averageIndex);
end

function change_plot(src,~)
    overlaySpectraOnMRI(src.String{src.Value});
end

function set_min(src, ~)
    max_obj = findobj('Tag', 'max');
    max_obj.String;
    ppm_min = str2double(src.String);
    ppm_max = str2double(max_obj.String);
    if(~isempty(ppm_min) && ~isempty(ppm_max))
        if(ppm_min <= ppm_max)
            plot_obj = findobj('Tag', 'plot_type');
            cur_type = plot_obj.String{plot_obj.Value};
            overlaySpectraOnMRI(cur_type, ppm_min, ppm_max);
        end
    end
end

function set_max(src, ~)
    min_obj = findobj('Tag', 'min');
    min_obj.String;
    ppm_max = str2double(src.String);
    ppm_min = str2double(min_obj.String);
    if(~isempty(ppm_min) && ~isempty(ppm_max))
        if(ppm_min <= ppm_max)
            plot_obj = findobj('Tag', 'plot_type');
            cur_type = plot_obj.String{plot_obj.Value};
            overlaySpectraOnMRI(cur_type, ppm_min, ppm_max);
        end
    end
end

function setMouseClickOption(src, ~)
    global st

    %get the axis object for the axial plane
    axis = st.vols{1}.ax{1}.ax;

    if(src.Value == 2)
        %If the popup menu is set to select_spec
        %set the callback function of click on axis to the method
        %select_spec.m
        set(axis, 'ButtonDownFcn', @select_spec);
        %position new axis here with properties
        axes('Visible','off', ...
             'Parent', st.fig, ...
             'YDir','normal', ...
             'Units','Pixels', ...
             'Box','on', ...
             'Position', [380 445 250 220], ...
             'Units','normalized', ...
             'Visible','on', ...
             'Tag', 'csi_axis');
    else
        %If the popup menu is set to reposition
        %set click on axial plane to reposition
        set(axis, 'ButtonDownFcn', @repos_start);
        %find the axis used to plot MRSI specs
        ax = findobj(gcf, 'Tag', 'csi_axis');
        %delete the specs
        delete(ax);
    end
end

%callback method used for repositioning
function repos_start(~, ~)
    currentFigure = gcbf;
    %checking if the click is a right click
    if ~strcmpi(get(currentFigure, 'SelectionType'),'alt')
        %set callback functions for movement of mouse and finished
        %clicking. spm_orthviews('reposition') is for movement of mouse
        % and repos_end for end of click.
        set(currentFigure, 'windowbuttonmotionfcn',@()(spm_orthviews('reposition')), 'windowbuttonupfcn',@repos_end);
        %call reposition in spm_orthviews
        spm_orthviews('reposition');
    end
end

function repos_end(~, ~)
    %remove callback functions for movement of mouse and wiindows
    %button up.
    set(gcbf,'windowbuttonmotionfcn','', 'windowbuttonupfcn','');
end

function showVoxelCallback(src, ~)
    global show_vox
    show_vox = src.Value;
    overlaySpectraOnMRI();
end

function checkArguments(in)
    if ~exist('spm_image', 'file')
        error('please add spm12 to your path before continuing');
    end
    checkSpatialFT(in)
end

function buildMouseClickDropdown(spmFigure)
    %Labels for dropdown
    labels={'Reposition', 'Select Spec'};
    %putting the dropdown menu in the figure to reposition or select spec
    mouseDropdown = uicontrol('Parent', spmFigure, ...
        'Units', 'Pixels', ...
        'Position', [320 400 100 20], ...
        'Style', 'Popupmenu', ...
        'String', labels);
    %selector callback function
    mouseDropdown.Callback = @setMouseClickOption;
end

function [maxPPMTextBox, minPPMTextBox] = buildPPMBoundTextBox(spmFigure)
    maxPPMTextBox = uicontrol('Parent', spmFigure, ...
        'Units', 'Pixels', ...
        'Position', [100, 405, 50, 20],...
        'Style', 'edit', ...
        'String', 'ppmmax', ...
        'Tag', 'max', ...
        'Callback', @set_max);

    minPPMTextBox = uicontrol('Parent', spmFigure, ...
        'Units', 'Pixels', ...
        'Position', [50, 405, 50, 20],...
        'Style', 'edit', ...
        'String', 'ppmmin', ...
        'Tag', 'min', ...
        'Callback', @set_min);
end

%
function buildPlotTypeDropdown(spmFigure)
    %labels for the dropdown menu
    plotLabels = {'real', 'imaginary', 'magnitude'};
    %creates a dropdown menu in the figure to decide to plot real, imaginary,
    %or magnitud
    plot_type_selector = uicontrol('Parent', spmFigure, ...
        'Units', 'Pixels', ...
        'Position', [200 400 100 20], ...
        'Style', 'Popupmenu', ...
        'String', plotLabels, ...
        'Tag', 'plot_type', ...
        'Callback', @change_plot);
end

function buildShowVoxelButton(spmFigure)
    uicontrol('Parent', spmFigure, ...
        'Units', 'Pixels', ...
        'Position', [500 405 60 20],...
        'Style', 'togglebutton', ...
        'String', 'Show Voxels', ...
        'Tag', 'show_voxels', ...
        'Callback', @showVoxelCallback);
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function varargout = OregonatorGUI(varargin)
% OREGONATORGUI M-file for OregonatorGUI.fig
%      OREGONATORGUI, by itself, creates a new OREGONATORGUI or raises the existing
%      singleton*.
%
%      H = OREGONATORGUI returns the handle to a new OREGONATORGUI or the handle to
%      the existing singleton*.
%
%      OREGONATORGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in OREGONATORGUI.M with the given input arguments.
%
%      OREGONATORGUI('Property','Value',...) creates a new OREGONATORGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before OregonatorGUI_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to OregonatorGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Copyright 2002-2003 The MathWorks, Inc.

% Edit the above text to modify the response to help OregonatorGUI

% Last Modified by GUIDE v2.5 02-Dec-2005 23:24:01

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @OregonatorGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @OregonatorGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before OregonatorGUI is made visible.
function OregonatorGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to OregonatorGUI (see VARARGIN)

% Choose default command line output for OregonatorGUI
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes OregonatorGUI wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = OregonatorGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

x0=[1 1 1];

[t,x]=ode15s(@oregonator,[0 600],x0);

axes(handles.axes1);
plot(t,x(:,1),'r');
axis([300,600,0,100]);
title('figure1');
Xlabel('temps');
Ylabel('x(1)');

axes(handles.axes2);
plot(t,x(:,2),'g');
axis([300,600,0,3]);
title('figure2');
Xlabel('temps');
Ylabel('x(2)');

axes(handles.axes3);
plot(t,x(:,3),'m');
axis([300,600,0,100]);
title('figure3');
Xlabel('temps');
Ylabel('x(3)');


text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu, J., Zheng, Y., Gao, Z., Jiang, Y., Hu, X., Zhu, Y., Jiao, J., Liu, M. (2020)
%           Quadratic Pose Estimation Problems: Unified Solutions,
%           Solvability/Observability Analysis and Uncertainty Description
%           in A Globally Optimal Framework.
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com

function [XX, ff] = optimize_quat_cov2(q, F, cov_left, solver, verbose)
if(~strcmp(solver, 'scs'))
    X = sdpvar(4, 4);
    f = q.' * X * q;
    cons = [
        X >= 0, ...
        sym2vec3(F * X * F.') == sym2vec3(cov_left)
        ];
    options = sdpsettings('solver', solver, 'verbose', verbose);
    if(strcmp(solver, 'sdpa_gmp'))
        options.sdpa_gmp.epsilonDash = 1.0e-35;
        options.sdpa_gmp.precision = 250;
    end
    optimize(cons, f, options);
    XX = value(X);
    ff = value(f);
else
    cvx_solver scs
    cvx_precision high

    cvx_begin sdp
         variable X(4,4) symmetric semidefinite
         dual variable Q
         minimize(q.' * X * q)
         X >= 0 * eye(4) : Q
         subject to
             sym2vec3(F * X * F.') == sym2vec3(cov_left)
    cvx_end
    XX = real(X);
    ff = 0;
end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% NANMIN Minimum value, ignoring NaNs.
%    M = NANMIN(A) returns the minimum of A with NaNs treated as missing.
%    For vectors, M is the smallest non-NaN element in A.  For matrices, M
%    is a row vector containing the minimum non-NaN element from each
%    column.  For N-D arrays, NANMIN operates along the first non-singleton
%    dimension.
%
%    [M,NDX] = NANMIN(A) returns the indices of the minimum values in A.  If
%    the values along the first non-singleton dimension contain more than
%    one minimal element, the index of the first one is returned.
%
%    M = NANMIN(A,B) returns an array the same size as A and B with the
%    smallest elements taken from A or B.  Either one can be a scalar.
%
%    [M,NDX] = NANMIN(A,[],DIM) operates along the dimension DIM.
%
%    See also MIN, NANMAX, NANMEAN, NANMEDIAN, NANVAR, NANSTD.
%
%    Reference page in Doc Center
%       doc stats/nanmin
%
%    Other functions named nanmin
%
%       distributed/nanmin    fints/nanmin
%

author: jmaih


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcRotate(Temp,Event,handles)
handles=guidata(handles.MU_matrix_display);

if ~isempty(handles.V.Segs)
    choice = questdlg('Segmentation mask is detected, transform operation will reset them, preceed?','Mask Reset', ...
                      'No, go save mask','Yes','No, go save mask');
    if isempty(choice)
        warndlg('Transform is cancelled.');
        return;
    end
    % Handle response
    switch choice
        case 'No, go save mask'
            warndlg('Save your mask before transform.');
            return;
    end

    handles.Mask=handles.Mask*0;
    handles.V.Segs=[];
end
% close 3D slicer
global Figure_handles
if isfield(Figure_handles,'MU_display2')
    slicer_display_handles=guidata(Figure_handles.MU_display2);
    if Figure_handles.MU_display == slicer_display_handles.Parent
        close(Figure_handles.MU_display2);
    end
end

if numel(handles.V.DimSize) == 3
    defaultOrigin = [ceil(handles.V.DimSize(2)/2) ceil(handles.V.DimSize(1)/2) ceil(handles.V.DimSize(3)/2)];
else
    defaultOrigin = [ceil(handles.V.DimSize(2)/2) ceil(handles.V.DimSize(1)/2) 0];
end

Input = inputdlg({'Please input a rotation origin [x,y,z].'; 'Please input a rotation axis [x,y,z].'; 'Please input a rotation angle (rad).'},'Input Value',1, ...
                 {['[' num2str(defaultOrigin(1)) ',' num2str(defaultOrigin(2)) ',' num2str(defaultOrigin(3)) ']'],'[0,0,1]','0'});
if isempty(Input)
    warndlg('Matrix rotation was cancelled.');
    return;
end

[Type,ok] = listdlg('ListString',{'Nearest','Linear','Cubic'}, ...
                    'SelectionMode','single',...
                    'PromptString','Interpolation Method',...
                    'Name','Interpolation');
if ok==0
    warndlg('Matrix rotation was cancelled.');
    return;
end

switch Type
    case 1
        interp = 'nearest';
    case 2
        interp = 'linear';
    case 3
        interp = 'cubic';
end

try
    MU_update_waitbar(handles.Progress_axes,1,3);
    pause(0.1);
    eval(['Ori=' Input{1} ';']);
    eval(['Axis=' Input{2} ';']);
    eval(['Angle=' Input{3} ';']);
    if Angle==0
        MU_update_waitbar(handles.Progress_axes,3,3);
        return;
    end
    [T, R]=rotate3DT_MU(Ori, Axis, Angle,interp);
    handles.TMatrix =rotate3D(handles.TMatrix, T, R);
    MU_update_waitbar(handles.Progress_axes,2,3);
    handles.Mask =round(rotate3D(handles.Mask, T, R));
    MergeM=get(handles.Matrix_name_edit,'String');
    set(handles.Matrix_name_edit,'String',[MergeM '_rot']);
    handles=MU_update_image(handles.Matrix_display_axes,{handles.TMatrix,handles.Mask},handles,0);
    MU_update_waitbar(handles.Progress_axes,3,3);
catch me
    errordlg('The input rotation value is invalid.');
    return;
end

guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
filename='BridgeCool_Quadrilateral_Bilinear_Structured';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = 1;
constraint = {'volume'};
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';

nsteps = 15;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%[2016]-"The whale optimization algorithm"

% (9/12/2020)

function WOA = jWhaleOptimizationAlgorithm(feat,label,opts)
% Parameters
lb    = 0;
ub    = 1;
thres = 0.5;
b     = 1;     % constant

if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'b'), b = opts.b; end
if isfield(opts,'thres'), thres = opts.thres; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = zeros(N,dim);
for i = 1:N
	for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
	end
end
% Fitness
fit  = zeros(1,N);
fitG = inf;
for i = 1:N
  fit(i) = fun(feat,label,(X(i,:) > thres),opts);
  % Global best
  if fit(i) < fitG
    fitG = fit(i);
    Xgb  = X(i,:);
  end
end
% Pre
curve = zeros(1,max_Iter);
curve(1) = fitG;
t = 2;
while t <= max_Iter
	% Define a, linearly decreases from 2 to 0
  a = 2 - t * (2 / max_Iter);
  for i = 1:N
    % Parameter A (2.3)
    A = 2 * a * rand() - a;
    % Paramater C (2.4)
    C = 2 * rand();
    % Parameter p, random number in [0,1]
    p = rand();
    % Parameter l, random number in [-1,1]
    l = -1 + 2 * rand();
    % Whale position update (2.6)
    if p  < 0.5
      % {1} Encircling prey
      if abs(A) < 1
        for d = 1:dim
          % Compute D (2.1)
          Dx     = abs(C * Xgb(d) - X(i,d));
          % Position update (2.2)
          X(i,d) = Xgb(d) - A * Dx;
        end
      % {2} Search for prey
      elseif abs(A) >= 1
        for d = 1:dim
          % Select a random whale
          k      = randi([1,N]);
          % Compute D (2.7)
          Dx     = abs(C * X(k,d) - X(i,d));
          % Position update (2.8)
          X(i,d) = X(k,d) - A * Dx;
        end
      end
    % {3} Bubble-net attacking
    elseif p >= 0.5
      for d = 1:dim
        % Distance of whale to prey
        dist   = abs(Xgb(d) - X(i,d));
        % Position update (2.5)
        X(i,d) = dist * exp(b * l) * cos(2 * pi * l) + Xgb(d);
      end
    end
    % Boundary
    XB = X(i,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
    X(i,:) = XB;
  end
  % Fitness
  for i = 1:N
    % Fitness
    fit(i) = fun(feat,label,(X(i,:) > thres),opts);
    % Global best
    if fit(i) < fitG
      fitG = fit(i);
      Xgb  = X(i,:);
    end
  end
  curve(t) = fitG;
  fprintf('\nIteration %d Best (WOA)= %f',t,curve(t))
  t = t + 1;
end
% Select features
Pos   = 1:dim;
Sf    = Pos((Xgb > thres) == 1);
sFeat = feat(:,Sf);
% Store results
WOA.sf = Sf;
WOA.ff = sFeat;
WOA.nf = length(Sf);
WOA.c  = curve;
WOA.f  = feat;
WOA.l  = label;
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
% std_getindvar - get independent variables of a STUDY
%
% Usage:
%   [indvar indvarvals] = std_getindvar(STUDY);
%   [indvar indvarvals] = std_getindvar(STUDY, mode, scandesign);
%
% Input:
%   STUDY - EEGLAB STUDY structure
%   mode  - ['datinfo'|'trialinfo'|'both'] get independent variables
%           linked to STUDY.datasetinfo, STUDY.datasetinfo.trialinfo or
%           both. Default is 'both'.
%   scandesign - [0|1] scan STUDY design for additional combinations of
%           independent variable values. Default is 0.
%
% Output:
%   indvar     - [cell array] cell array of independent variable names
%   indvarvals - [cell array] cell array of independent variable values
%
% Authors: A. Delorme, CERCO/CNRS and SCCN/UCSD, 2010-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, 2010, arno@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [ factor, factorvals, subjects, paired ] = std_getindvar(STUDY, mode, scandesign)

if nargin < 1
    help std_getindvar;
    return;
end

factor     = {};
factorvals = {};
paired     = {};
subjects   = {};
if nargin < 2, mode = 'both'; end
if nargin < 3, scandesign = 0; end
countfact = 1;
setinfo = STUDY.datasetinfo;

if strcmpi(mode, 'datinfo') || strcmpi(mode, 'both')
    % get trial info
    ff = fieldnames(setinfo);
    ff = setdiff_bc(ff, { 'filepath' 'filename' 'subject' 'index' 'ncomps' 'comps' 'trialinfo' });
    for index = 1:length(ff)
        % check we have the same data type
        allSetinfoVals = { setinfo.(ff{index}) };
        if any(cellfun(@ischar, allSetinfoVals))
            if ~all(cellfun(@ischar, allSetinfoVals))
                allSetinfoVals = cellfun(@num2str, allSetinfoVals, 'uniformoutput', false);
            end
        end

        if ischar(allSetinfoVals{1})
            tmpvals = unique_bc(allSetinfoVals);
            if length(tmpvals) > 1
                factor{    countfact} = ff{index};
                factorvals{countfact} = tmpvals;

                % get subject for each factor value
                intersectSubject = { setinfo(:).subject };
                for c = 1:length(tmpvals)
                    datind = strmatch(tmpvals{c}, allSetinfoVals, 'exact');
                    subjects{  countfact}{c} = unique_bc( { setinfo(datind).subject } );
                    intersectSubject = intersect(intersectSubject, subjects{  countfact}{c});
                end

                numValues = cellfun(@length, subjects{  countfact});
                if length(intersectSubject) == numValues(1) && length(unique(numValues)) == 1
                    paired{countfact} = 'on'; % trial data always paired
                else
                    paired{countfact} = 'off';
                end
                countfact = countfact + 1;
            end
        else
            tmpvals = unique_bc([ allSetinfoVals{:} ] );
            if length(tmpvals) > 1
                factor{    countfact} = ff{index};
                factorvals{countfact} = mattocell(tmpvals);

                % get subject for each factor value
                intersectSubject = { setinfo(:).subject };
                for c = 1:length(tmpvals)
                    eval( [ 'datind = find(tmpvals(c) == [ setinfo.' ff{index} ']);' ] );
                    subjects{  countfact}{c} = unique_bc( { setinfo(datind).subject } );
                    intersectSubject = intersect(intersectSubject, subjects{  countfact}{c});
                end

                numValues = cellfun(@length, subjects{  countfact});
                if length(intersectSubject) == numValues(1) && length(unique(numValues)) == 1
                    paired{countfact} = 'on'; % trial data always paired
                else
                    paired{countfact} = 'off';
                end
                countfact = countfact + 1;
            end
        end
    end
end

% add ind. variables for trials
% -----------------------------
if strcmpi(mode, 'trialinfo') || strcmpi(mode, 'both')
    % add trial info
    if isfield(setinfo, 'trialinfo')
        ff = fieldnames(setinfo(1).trialinfo);
        for index = 1:length(ff)
            % check if any of the datasets are using string for event type
            allFieldsPresent = cellfun(@(x)(isfield(x, ff{index})), { setinfo.trialinfo });
            allFirstVal = cellfun(@(x)(getfield(x, ff{index})), { setinfo(allFieldsPresent).trialinfo }, 'uniformoutput', false);

            if any(cellfun(@isstr, allFirstVal))
                alltmpvals = {};
                for ind = 1:length(setinfo)
                    if isfield(setinfo(ind).trialinfo, ff{index})
                        eval( [ 'tmpTrialVals = { setinfo(ind).trialinfo.' ff{index} ' };' ] );
                        if isnumeric(tmpTrialVals{1})
                            % convert to string if necessary
                            tmpTrialVals = cellfun(@num2str, tmpTrialVals, 'uniformoutput', false);
                        end
                        tmpvals = unique_bc(tmpTrialVals);
                    else tmpvals = {};
                    end
                    if isempty(alltmpvals)
                        alltmpvals = tmpvals;
                    else
                        alltmpvals = { alltmpvals{:} tmpvals{:} };
                    end
                end
                alltmpvals = unique_bc(alltmpvals);
                if length(alltmpvals) > 1
                    factor{    countfact} = ff{index};
                    factorvals{countfact} = alltmpvals;
                    subjects{  countfact} = {};
                    paired{countfact}     = 'on'; % trial data always paired
                    countfact = countfact + 1;
                end
            else
                alltmpvals = [];
                for ind = 1:length(setinfo)
                    if isfield(setinfo(ind).trialinfo, ff{index}) && ~iscell( setinfo(ind).trialinfo(1).(ff{index}))
                        eval( [ 'tmpvals = unique_bc([ setinfo(ind).trialinfo.' ff{index} ' ]);' ] );
                    else tmpvals = [];
                    end
                    alltmpvals = [ alltmpvals tmpvals ];
                end
                alltmpvals = unique_bc(alltmpvals);
                if length(alltmpvals) > 1
                    factor{    countfact} = ff{index};
                    factorvals{countfact} = mattocell(alltmpvals);
                    subjects{  countfact} = {};
                    paired{countfact}     = 'off'; % pairing is irrelevant for continuous var
                    countfact = countfact + 1;
                end
            end
        end
    end
end

% scan existing design for additional combinations
% ------------------------------------------------
if scandesign
    for desind = 1:length(STUDY.design)
        for iVar = 1:length(STUDY.design(desind).variable)
            if strcmpi(STUDY.design(desind).variable(iVar).vartype, 'categorical')
                pos1 = strmatch(STUDY.design(desind).variable(iVar).label, factor, 'exact');
                if ~isempty(pos1), add1 = mysetdiff(STUDY.design(desind).variable(iVar).value, factorvals{pos1}); else add1 = []; end
                for addInd = 1:length(add1)
                    duplicate = 0;
                    for iVarVal = 1:length(factorvals{pos1})
                        if isequal(factorvals{pos1}{iVarVal}, add1{addInd})
                            duplicate = 1;
                        end
                    end
                    if ~duplicate
                        factorvals{pos1} = { factorvals{pos1}{:} add1{addInd} };
                    end
                end
            end
        end
    end
end

function cellout = mysetdiff(cell1, cell2);

    if ischar(cell2{1})
         indcell = cellfun(@iscell, cell1);
    else indcell = cellfun(@(x)(length(x)>1), cell1);
    end
    cellout = cell1(indcell);

%     if ischar(cell1{1}) && ischar(cell2{1})
%          cellout = setdiff_bc(cell1, cell2);
%     elseif ~ischar(cell1{1}) && ~ischar(cell2{1})
%         cellout = mattocell(setdiff( [ cell1{:} ], [ cell2{:} ]));
%     elseif ischar(cell1{1}) && ~ischar(cell2{1})
%          cellout = setdiff_bc(cell1, cellfun(@(x)(num2str(x)),cell2, 'uniformoutput', false));
%     elseif ~ischar(cell1{1}) && ischar(cell2{1})
%          cellout = setdiff_bc(cellfun(@(x)(num2str(x)),cell1, 'uniformoutput', false), cell2);
%     end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function train_base_net_x4(varargin)
% load data & label
data = load('./data/train/SDR_youtube_80_x4.mat') ;
label = load('./data/train/HDR_youtube_80.mat') ;
imdb.images.data = data.SDR_data;
imdb.images.label = label.HDR_data;
imdb.images.set = cat(2, ones(1, size(data.SDR_data, 4)-500), 2*ones(1, 500));

% set CNN model
net = net_base_x4();

% set the learning rate and weight decay for biases
% default values are used for filters
for i = 2:2:114
    net.params(i).learningRate = 0.1;
    net.params(i).weightDecay = 0;
end
net.conserveMemory = true;

% options
opts.solver = @adam;
opts.train.batchSize = 16;
opts.train.continue = true;
opts.train.gpus = 1;
opts.train.prefetch = false ;
opts.train.expDir = './net/net_base_x4' ;
opts.train.learningRate = 5*1e-7*ones(1, 200);
opts.train.weightDecay = 0.0005;
opts.train.numEpochs = numel(opts.train.learningRate) ;
opts.train.derOutputs = {'objective', 1} ;
[opts, ~] = vl_argparse(opts.train, varargin) ;

% record
if(~isdir(opts.expDir))
    mkdir(opts.expDir);
end

% call training function
[net, info] = cnn_train_dag(net, imdb, @getBatch, opts) ;

function inputs = getBatch(imdb, batch, opts)
image = imdb.images.data(:, :, :, batch) ;
label = imdb.images.label(:, :, :, batch) ;

image = single(image)/255;
label = single(label)/1023;
inputs = {'input', gpuArray(image), 'label', gpuArray(label)};

author: sooyekim

repo: Deep-SR-ITM


text:
% BS2RV.m - Binary string to real vector
%
% This function decodes binary chromosomes into vectors of reals. The
% chromosomes are seen as the concatenation of binary strings of given
% length, and decoded into real numbers in a specified interval using
% either standard binary or Gray decoding.
%
% Syntax:       Phen = bs2rv(Chrom,FieldD)
%
% Input parameters:
%
%               Chrom    - Matrix containing the chromosomes of the current
%                          population. Each line corresponds to one
%                          individual's concatenated binary string
%			   representation. Leftmost bits are MSb and
%			   rightmost are LSb.
%
%               FieldD   - Matrix describing the length and how to decode
%			   each substring in the chromosome. It has the
%			   following structure:
%
%				[len;		(num)
%				 lb;		(num)
%				 ub;		(num)
%				 code;		(0=binary     | 1=gray)
%				 scale;		(0=arithmetic | 1=logarithmic)
%				 lbin;		(0=excluded   | 1=included)
%				 ubin];		(0=excluded   | 1=included)
%
%			   where
%				len   - row vector containing the length of
%					each substring in Chrom. sum(len)
%					should equal the individual length.
%				lb,
%				ub    - Lower and upper bounds for each
%					variable.
%				code  - binary row vector indicating how each
%					substring is to be decoded.
%				scale - binary row vector indicating where to
%					use arithmetic and/or logarithmic
%					scaling.
%				lbin,
%				ubin  - binary row vectors indicating whether
%					or not to include each bound in the
%					representation range
%
% Output parameter:
%
%               Phen     - Real matrix containing the population phenotypes.

%
% Author: Carlos Fonseca, 	Updated: Andrew Chipperfield
% Date: 08/06/93,		Date: 26-Jan-94

function Phen = bs2rv(Chrom,FieldD)

% Identify the population size (Nind)
%      and the chromosome length (Lind)
[Nind,Lind] = size(Chrom);

% Identify the number of decision variables (Nvar)
[seven,Nvar] = size(FieldD);

if seven ~= 7
	error('FieldD must have 7 rows.');
end

% Get substring properties
len = FieldD(1,:);
lb = FieldD(2,:);
ub = FieldD(3,:);
code = ~(~FieldD(4,:));
scale = ~(~FieldD(5,:));
lin = ~(~FieldD(6,:));
uin = ~(~FieldD(7,:));

% Check substring properties for consistency
if sum(len) ~= Lind,
	error('Data in FieldD must agree with chromosome length');
end

if ~all(lb(scale).*ub(scale)>0)
	error('Log-scaled variables must not include 0 in their range');
end

% Decode chromosomes
Phen = zeros(Nind,Nvar);

lf = cumsum(len);
li = cumsum([1 len]);
Prec = .5 .^ len;

logsgn = sign(lb(scale));
lb(scale) = log( abs(lb(scale)) );
ub(scale) = log( abs(ub(scale)) );
delta = ub - lb;

Prec = .5 .^ len;
num = (~lin) .* Prec;
den = (lin + uin - 1) .* Prec;

for i = 1:Nvar,
    idx = li(i):lf(i);
    if code(i) % Gray decoding
	    Chrom(:,idx)=rem(cumsum(Chrom(:,idx)')',2);
    end
    Phen(:,i) = Chrom(:,idx) * [ (.5).^(1:len(i))' ];
    Phen(:,i) = lb(i) + delta(i) * (Phen(:,i) + num(i)) ./ (1 - den(i));
end

expand = ones(Nind,1);
if any(scale)
	Phen(:,scale) = logsgn(expand,:) .* exp(Phen(:,scale));
end

author: vonsylvia


text:
function VLMplay(varargin)
% Visualizes the polygons from a video and the tracks if available inside the annotation.
%
% To play the annotations only:
% VLMplay(annotation)
%
% To get video for selected frames
% VLMplay(D, ndx, HOMEVIDEOS, frameSkipRate)

switch length(varargin)
    case 1
        v = varargin{1};
        annotation = v.annotation;
        videoFramesPath = '';

    case 3
        D = varargin{1};
        ndx = varargin{2};

        if ndx>length(D)
            error('Index outside database range')
        end

        HOMEVIDEOS = varargin{3};

        annotation = D(ndx).annotation;
        videoFramesPath = fullfile(HOMEVIDEOS, annotation.folder, annotation.filename(1:end-4));

    case 4
        D = varargin{1};
        ndx = varargin{2};
        if ndx>length(D)
            error('Index outside database range')
        end

        HOMEVIDEOS = varargin{3};

        annotation = D(ndx).annotation;
        videoFramesPath = fullfile(HOMEVIDEOS, annotation.folder, annotation.filename(1:end-4));
        frameSkipRate = varargin{4};

    otherwise
        error('wrong number of input arguments.')
end

numFrames = str2num(annotation.numFrames);
% Define colors
colors = hsv(15);

if ~isempty(videoFramesPath)
    % Open the directory containing the frames
    videoDir = dir(fullfile(videoFramesPath, '*.jpg'));
end

%default if frameIdxs is not given
if(~exist('frameSkipRate', 'var'))
    frameSkipRate = 1;
end

% Draw image
figure(1);
% if we want to resize, uncomment the next line:
%set(gcf, 'position', [183   712   202   167]);
rootName = annotation.filename(1:end-4);

for f = 1:frameSkipRate:numFrames
    [annotation, frame] = VLMgetvideoframe(D(ndx).annotation, f, HOMEVIDEOS);
    [nrows ncols cc] = size(frame);
    hold off;
    LMplot(annotation, frame)

    hold on
    plot([annotation.tracks.track.x], [annotation.tracks.track.y], 'o', 'MarkerFaceColor','r', 'MarkerEdgeColor','y', 'MarkerSize', 5)

    axis([1 ncols 1 nrows])
    drawnow
end

author: CSAILVision

repo: LabelMeToolbox


text:
% RPC processes requests to the slaves.  It is the main mechanism by which
% the server communicates with the slaves and has a "remote procedure
% call"-like interface.
%
% rpc(reqHook, ...) runs the function identified by the string 'reqHook' on
% all of the slaves with the arguments that follow.
%
% The arguments that follow reqHook are passed in the same order as the
% arguments would be passed to the function if it were called locally.
% However, each argument represents an array of arguments passed to the
% various slaves.  In particular, for a server with N slaves, each argument
% is expected to be a cell array of N entries, where the i'th entry of the
% array contains the argument given to the i'th slave.  Several additional
% special cases are handled for convienience: (i) If the argument is not a
% cell array, it is wrapped in a cell array by arg={arg}, and (ii) If the
% argument is a singleton (as produced by (i), or is a cell array with only
% 1 element), the argument is replicated for all N slaves.
%
% Thus, for a 2 slave system, the following function calls are all identical:
%
% server.rpc('zeros', 3,1);
% server.rpc('zeros', {3},{1});
% server.rpc('zeros', {3,3},{1,1});
%
% [result1, ...] = server.rpc(...) works as above, but returns in
% result1, ... the results of the remote function calls.  The results are
% in the same order as the results returned by the remote function, but as
% with arguments, they are cell arrays containing a single result per slave.
%
% Example for 2 slaves:
%
% Z = server.rpc('zeros', 1,3);
% Z
%    Z = [ 1x3 double ]   [ 1x3 double ]
% Z{1}
%    ans = 0  0  0
%
% Caveat:
%
% At present, if server.rpc is invoked without output arguments
% the system will not retrieve any outputs from the remote function call.
% Thus, whereas MATLAB normally returns the first output (and prints the
% result) when no outputs are given this function will return nothing.
%
function varargout = rpc(server, reqHook, varargin)

  if (nargin >= 3)

    % check arguments and convert type if necessary
    for j=1:length(varargin)
      if (length(varargin{j}) == 0)
        error(['Empty argument array for argument ' num2str(j)]);
      end
      if (~isa(varargin{j}, 'cell'))
        varargin{j} = {varargin{j}};
      end
    end

    % rearrange arguments to be per-slave
    for i=1:server.slaveCount
      for j=1:length(varargin)
        if (length(varargin{j}) == 1)
          slaveVars{i}{j} = varargin{j}{1};
        else
          slaveVars{i}{j} = varargin{j}{i};
        end
      end
    end

    % send request to each slave
    for i=1:server.slaveCount
      if (~isempty(slaveVars{i}))
        server.sendRequest(i, nargout, reqHook, slaveVars{i}{:});
      else
        server.sendRequest(i, nargout, reqHook);
      end
    end

  else

    % send request to each slave;  no args
    for i=1:server.slaveCount
      server.sendRequest(i, nargout, reqHook);
    end
  end

  % get replies from each slave
  errCount=0;
  for i=1:server.slaveCount
    try
      slaveResults{i}=server.getReply(i);
    catch
      err=lasterror;
      err.message = sprintf('Error from slave %d:\n%s', i, err.message);
      errCount = errCount+1;
    end
  end
  if (errCount > 0)
    err.message = sprintf('%s\n\n %d slaves returned errors.', err.message, errCount);
    rethrow(err);
  end

  % reformat result
  for j=1:nargout
    varargout{j}={};
    for i=1:server.slaveCount
      varargout{j}{i} = slaveResults{i}{j};
    end
  end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function [ a, ipvt, rcond ] = chico ( a, lda, n )

%*****************************************************************************80
%
%% CHICO factors a complex hermitian matrix and estimates its condition.
%
%  Discussion:
%
%    If RCOND is not needed, CHIFA is slightly faster.
%
%    To solve A*X = B, follow CHICO by CHISL.
%
%    To compute inverse(A)*C, follow CHICO by CHISL.
%
%    To compute inverse(A), follow CHICO by CHIDI.
%
%    To compute determinant(A), follow CHICO by CHIDI.
%
%    To compute inertia(A), follow CHICO by CHIDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 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 A(LDA,N); on input, the hermitian matrix to be
%    factored.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex A(LDA,N); a block diagonal matrix and the multipliers
%    which were used to obtain it.  The factorization can be written
%    A = U*D*hermitian(U) where U is a product of permutation and unit
%    upper triangular matrices, hermitian(U) is the conjugate transpose
%    of U, and D is block diagonal with 1 by 1 and 2 by 2 blocks.
%    Only the diagonal and upper triangle are used.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, real RCOND, an estimate of the reciprocal condition of
%    the matrix.  For the system A*X = B, relative perturbations in A and B
%    of size EPSILON may cause relative perturbations in X of size
%    (EPSILON/RCOND).  If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate underflows.
%
%  Local Parameters:
%
%    Local, complex Z(N), a work vector whose contents are usually
%    unimportant.  If A is close to a singular matrix, then Z is an
%    approximate null vector in the sense that
%      norm(A*Z) = RCOND * norm(A) * norm(Z).
%

%
%  Find norm of A using only upper half.
%
  for j = 1 : n

    z(j) = scasum ( j, a(1:j,j), 1 );

    for i = 1 : j - 1
      z(i) = real ( z(i) ) + cabs1 ( a(i,j) );
    end

  end

  anorm = 0.0;
  for j = 1 : n
    anorm = max ( anorm, real ( z(j) ) );
  end
%
%  Factor.
%
  [ a, ipvt, info ] = chifa ( a, lda, n );
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A*Z = Y and A*Y = E.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where U*D*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve U*D*W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  k = n;

  while ( 0 < k )

    if ( ipvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    kp = abs ( ipvt(k) );
    kps = k + 1 - ks;

    if ( kp ~= kps )
      t      = z(kps);
      z(kps) = z(kp);
      z(kp)  = t;
    end

    if ( cabs1 ( z(k) ) ~= 0.0 )
      ek = csign1 ( ek, z(k) );
    end

    z(k) = z(k) + ek;
    z(1:k-ks) = z(1:k-ks) + z(k) * transpose ( a(1:k-ks,k) );

    if ( ks ~= 1 )

      if ( cabs1 ( z(k-1) ) ~= 0.0 )
        ek = csign1 ( ek, z(k-1) );
      end

      z(k-1) = z(k-1) + ek;
      z(1:k-ks) = z(1:k-ks) + z(k-1) * a(1:k-ks,k-1);

    end

    if ( ks ~= 2 )

      if ( cabs1 ( a(k,k) ) < cabs1 ( z(k) ) )
        s = cabs1 ( a(k,k) ) / cabs1 ( z(k) );
        z(1:n) = z(1:n) * s;
        ek = s * ek;
      end

      if ( cabs1 ( a(k,k) ) ~= 0.0 )
        z(k) = z(k) / a(k,k);
      else
        z(k) = 1.0;
      end

    else

      ak = a(k,k) / conj ( a(k-1,k) );
      akm1 = a(k-1,k-1) / a(k-1,k);
      bk = z(k) / conj ( a(k-1,k) );
      bkm1 = z(k-1) / a(k-1,k);
      denom = ak * akm1 - 1.0;
      z(k) = ( akm1 * bk - bkm1 ) / denom;
      z(k-1) = ( ak * bkm1 - bk ) / denom;

    end

    k = k - ks;

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
%
%  Solve hermitian(U) * Y = W.
%
  k = 1;

  while ( k <= n )

    if ( ipvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    if ( k ~= 1 )

      z(k) = z(k) + z(1:k-1) * conj ( a(1:k-1,k) );

      if ( ks == 2 )
        z(k+1) = z(k+1) + z(1:k-1) * conj ( a(1:k-1,k+1) );
      end

      kp = abs ( ipvt(k) );

      if ( kp ~= k )
        t     = z(k);
        z(k)  = z(kp);
        z(kp) = t;
      end

    end

    k = k + ks;

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = 1.0;
%
%  Solve U*D*V = Y.
%
  k = n;

  while ( 0 < k )

    if ( ipvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    if ( k ~= ks )

      kp = abs ( ipvt(k) );
      kps = k + 1 - ks;

      if ( kp ~= kps )
        t      = z(kps);
        z(kps) = z(kp);
        z(kp)  = t;
      end

      z(1:k-ks) = z(1:k-ks) + z(k) * transpose ( a(1:k-ks,k) );

      if ( ks == 2 )
        z(1:k-ks) = z(1:k-ks) + z(k-1) * a(1:k-ks,k-1);
      end

    end

    if ( ks ~= 2 )

      if ( cabs1 ( a(k,k) ) < cabs1 ( z(k) ) )
        s = cabs1 ( a(k,k) ) / cabs1 ( z(k) );
        z(1:n) = z(1:n) * s;
        ynorm = s * ynorm;
      end

      if ( cabs1 ( a(k,k) ) ~= 0.0 )
        z(k) = z(k) / a(k,k);
      else
        z(k) = 1.0;
      end

    else

      ak = a(k,k) / conj ( a(k-1,k) );
      akm1 = a(k-1,k-1) / a(k-1,k);
      bk = z(k) / conj ( a(k-1,k) );
      bkm1 = z(k-1) / a(k-1,k);
      denom = ak * akm1 - 1.0;
      z(k) = ( akm1 * bk - bkm1 ) / denom;
      z(k-1) = ( ak * bkm1 - bk ) / denom;

    end

    k = k - ks;

  end

  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;
%
%  Solve hermitian(U) * Z = V.
%
  k = 1;

  while ( k <= n )

    if ( ipvt(k) < 0 )
      ks = 2;
    else
      ks = 1;
    end

    if ( k ~= 1 )

      z(k) = z(k) + z(1:k-1) * conj ( a(1:k-1,k) );

      if ( ks == 2 )
        z(k+1) = z(k+1) + z(1:k-1) *  conj ( a(1:k-1,k+1) );
      end

      kp = abs ( ipvt(k) );

      if ( kp ~= k )
        t     = z(k);
        z(k)  = z(kp);
        z(kp) = t;
      end

    end

    k = k + ks;

  end
%
%  Make ZNORM = 1.
%
  s = 1.0 / scasum ( n, z, 1 );
  z(1:n) = z(1:n) * s;
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p = PlotHVLines(positions,direction,varargin)

%PlotHVLines - Plot vertical (resp. horizontal) lines at listed x (resp. y).
%
%  USAGE
%
%    p = PlotHVLines(positions,direction,options)
%
%    positions      list of abscissae/ordinates
%    direction      optional direction: 'h' or 'v' (default = 'v')
%    <options>      options for function plot
%

% Copyright (C) 2008-2012 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.

if nargin < 1,
 	error('Incorrect number of parameters (type ''help PlotHVLines'' for details).');
end
if nargin < 2,
	direction = 'v';
else
	direction = lower(direction);
end
if min(size(positions)) > 2,
 	error('List of abscissae/ordinates is not a vector (type ''help PlotHVLines'' for details).');
else
	positions = positions(:);
end

if ~isstring_FMAT(direction,'h','v'),
	varargin = {direction,varargin{:}};
	direction = 'v';
end

hold on;
if strcmp(direction,'v'),
	yLim = ylim;
	for i = 1:size(positions,1),
		plot([positions(i,1) positions(i,1)],yLim,varargin{:});
	end
else
	xLim = xlim;
	for i = 1:size(positions,1),
		plot(xLim,[positions(i,1) positions(i,1)],varargin{:});
	end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function g_nm = conjCoeffs(f_nm)
%CONJCOEFFS Get the complex SH coefficients of a conjugate sph. function
%
%   f_nm:   (N+1)^2 coefficients of original spherical function f(\theta,\phi)
%
%   g_nm:   (N+1)^2 coefficients of conjugate spherical function
%               g(\theta,\phi) = (f(\theta,\phi))*
%
%   The conversion is based on the relation Y^*_{nm} = (-1)^m Y_{n(-m)}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Archontis Politis, 10/10/2013
%   archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% order
N = sqrt(length(f_nm)) - 1;

g_nm = zeros((N+1)^2, 1);
for n=0:N
    for m=-n:n
        q_g = n*(n+1)+m;
        q_f = n*(n+1)-m;
        g_nm(q_g+1) = (-1)^m * conj(f_nm(q_f+1));
    end
end

end

author: polarch

repo: Spherical-Harmonic-Transform


text:
function x = gumbel_cdf_inv ( cdf )

%*****************************************************************************80
%
%% GUMBEL_CDF_INV inverts the Gumbel CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Output, real X, the corresponding argument of the CDF.
%
  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GUMBEL_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'GUMBEL_CDF_INV - Fatal error!' );
  end

  x =  - log ( - log ( cdf ) );

  return
end


text:
function trackers = config_tracker()
% config trackers to be evaluated
trackers = {struct('name', 'BACF',       'publish', 'ICCV-17') ...
            struct('name', 'CFNet',      'publish', 'CVPR-17') ...
            struct('name', 'CN',         'publish', 'CVPR-14') ...
            struct('name', 'CSK',        'publish', 'ECCV-12') ...
            struct('name', 'CSRDCF',     'publish', 'CVPR-17') ...
            struct('name', 'ASLA',       'publish', 'CVPR-12') ...
            struct('name', 'CT',         'publish', 'ECCV-12') ...
            struct('name', 'DSiam',      'publish', 'ICCV-17') ...
            struct('name', 'DSST',       'publish', 'BMVC-14') ...
            struct('name', 'ECO_HC',     'publish', 'CVPR-17') ...
            struct('name', 'ECO',        'publish', 'CVPR-17') ...
            struct('name', 'HCFT',       'publish', 'ICCV-15') ...
            struct('name', 'fDSST',      'publish', 'TPAMI-17')...
            struct('name', 'PTAV',       'publish', 'ICCV-17')...
            struct('name', 'IVT',        'publish', 'IJCV-08') ...
            struct('name', 'KCF',        'publish', 'TPAMI-15')...
            struct('name', 'L1APG',      'publish', 'CVPR-12') ...
            struct('name', 'LCT',        'publish', 'CVPR-15') ...
            struct('name', 'MEEM',       'publish', 'ECCV-14') ...
            struct('name', 'MIL',        'publish', 'CVPR-09') ...
            struct('name', 'SAMF',       'publish', 'ECCVW-14')...
            struct('name', 'SCT4',       'publish', 'CVPR-16') ...
            struct('name', 'SiamFC',     'publish', 'ECCVW-16')...
            struct('name', 'SINT',       'publish', 'CVPR-16') ...
            struct('name', 'SRDCF',      'publish', 'ICCV-15') ...
            struct('name', 'Staple_CA',  'publish', 'CVPR-17')...
            struct('name', 'Staple',     'publish', 'CVPR-16')...
            struct('name', 'STC',        'publish', 'ECCV-14') ...
            struct('name', 'STRCF',      'publish', 'CVPR-18') ...
            struct('name', 'Struck',     'publish', 'TPAMI-16')...
            struct('name', 'TLD',        'publish', 'TPAMI-12')...
            struct('name', 'MDNet',      'publish', 'CVPR-16')...
            struct('name', 'StructSiam', 'publish', 'ECCV-18')...
            struct('name', 'VITAL',      'publish', 'CVPR-18')...
            struct('name', 'TRACA',      'publish', 'CVPR-18')};
end

author: HengLan


text:
clear;

%%
% Find the dataset location
if exist('D:/Datasets/HeadPose', 'file')
    database_root = 'D:/Datasets/HeadPose/';
elseif(exist([getenv('USERPROFILE') '/Dropbox/AAM/test data/'], 'file'))
    database_root = [getenv('USERPROFILE') '/Dropbox/AAM/test data/'];
elseif(exist(['F:/Dropbox/Dropbox/AAM/test data/'], 'file'))
    database_root = 'F:/Dropbox/Dropbox/AAM/test data/';
elseif(exist(['F:/Dropbox/AAM/test data/'], 'file'))
    database_root = 'F:/Dropbox/AAM/test data/';
elseif(exist('/media/tadas/5E08AE0D08ADE3ED/Dropbox/AAM/test data', 'file'))
    database_root = '/media/tadas/5E08AE0D08ADE3ED/Dropbox/AAM/test data';
end

buDir = [database_root, '/bu/uniform-light/'];

% The fast and accurate clnf
%%
[resFolderBU_OF] = run_bu_experiment(buDir, false, 'model', 'model/main_clnf_general.txt');
[bu_error_OF, pred_hp_bu, gt_hp_bu, all_errors_bu_OF, rels_bu] = calcBUerror(resFolderBU_OF, buDir);

%%
% Run the Biwi test
biwi_dir = '/biwi pose/';

[res_folder_biwi_OF] = run_biwi_experiment(database_root, biwi_dir, false, false, 'model', 'model/main_clnf_general.txt');
% Calculate the resulting errors
[biwi_error_OF, pred_hp_biwi, gt_hp_biwi, ~, all_errors_biwi_OF, rels_biwi] = calcBiwiError(res_folder_biwi_OF, [database_root biwi_dir]);

%% Run the ICT test
ict_dir = ['/ict/'];

% Intensity
[res_folder_ict_OF] = run_ict_experiment(database_root, ict_dir, false, false, 'model', 'model/main_clnf_general.txt');
% Calculate the resulting errors
[ict_error_OF, pred_hp_ict, gt_hp_ict, ~, all_errors_ict_OF, rel_ict] = calcIctError(res_folder_ict_OF, [database_root ict_dir]);

%% Save the results
filename = 'results/Pose_OF_CLNF';
save(filename);

% Also save them in a reasonable .txt format for easy comparison
f = fopen('results/Pose_OF_CLNF.txt', 'w');
fprintf(f, 'Dataset and model,        pitch,  yaw,  roll,  mean,  median\n');
fprintf(f, 'biwi error:  %.3f,   %.3f, %.3f,  %.3f,  %.3f\n', biwi_error_OF, mean(all_errors_biwi_OF(:)), median(all_errors_biwi_OF(:)));
fprintf(f, 'bu error:    %.3f,   %.3f, %.3f,  %.3f,  %.3f\n', bu_error_OF, mean(all_errors_bu_OF(:)), median(all_errors_bu_OF(:)));
fprintf(f, 'ict error:   %.3f,   %.3f, %.3f,  %.3f,  %.3f\n', ict_error_OF, mean(all_errors_ict_OF(:)), median(all_errors_ict_OF(:)));

fclose(f);
clear 'f'

author: TadasBaltrusaitis

repo: OpenFace


text:
function [S2G, r]= refine(S2G)
% refine S2Grid
%
% Input
%  S2G - @S2Grid
%
% Output
%  S2G - @S2Grid with half the resolution

S2G.opt.resolution = S2G.opt.resolution / 2;
S2G.thetaGrid = refine(S2G.thetaGrid);
for i=1:length(S2G.rhoGrid)
  S2G.rhoGrid(2*i-1) = refine(S2G.rhoGrid(i));
end

theta = double(S2G.thetaGrid);
for i = 2:2:GridLength(S2G.thetaGrid)

  % construct full circle
  nrho= max(round(sin(theta(i)) * 2*pi/S2G.res),1);
  rho = (0:nrho-1)*2*pi/nrho;

  % take only those that are close to existing ones
  d = min(dist_outer(S2G.rhoGrid(i-1),rho));
  if i+1 <= length(theta)
    d = min(d,min(dist_outer(S2G.rhoGrid(i+1),rho)));
  end

  rho = rho(sin(theta(i-1))*d < S2G.res*1.5);
  S2G.rhoGrid(i) = S1Grid(rho,...
    getMin(S2G.rhoGrid(i-1)),getMax(S2G.rhoGrid(i-1)),check_option(S2G.rhoGrid(i-1)));

end

[S2G.x,S2G.y,S2G.z] = double(calcGrid(S2G.thetaGrid,S2G.rhoGrid));

end

author: mtex-toolbox

repo: mtex

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