text:
function vl_compilenn(varargin)
%VL_COMPILENN Compile the MatConvNet toolbox.
%   The `vl_compilenn()` function compiles the MEX files in the
%   MatConvNet toolbox. See below for the requirements for compiling
%   CPU and GPU code, respectively.
%
%   `vl_compilenn('OPTION', ARG, ...)` accepts the following options:
%
%   `EnableGpu`:: `false`
%      Set to true in order to enable GPU support.
%
%   `Verbose`:: 0
%      Set the verbosity level (0, 1 or 2).
%
%   `Debug`:: `false`
%      Set to true to compile the binaries with debugging
%      information.
%
%   `CudaMethod`:: Linux & Mac OS X: `mex`; Windows: `nvcc`
%      Choose the method used to compile the CUDA code. There are two
%      methods:
%
%      * The **`mex`** method uses the MATLAB MEX command with the
%        configuration file
%        `<MatConvNet>/matlab/src/config/mex_CUDA_<arch>.[sh/xml]`
%        This configuration file is in XML format since MATLAB 8.3
%        (R2014a) and is a Shell script for earlier versions. This
%        is, principle, the preferred method as it uses the
%        MATLAB-sanctioned compiler options.
%
%      * The **`nvcc`** method calls the NVIDIA CUDA compiler `nvcc`
%        directly to compile CUDA source code into object files.
%
%        This method allows to use a CUDA toolkit version that is not
%        the one that officially supported by a particular MATALB
%        version (see below). It is also the default method for
%        compilation under Windows and with CuDNN.
%
%   `CudaRoot`:: guessed automatically
%      This option specifies the path to the CUDA toolkit to use for
%      compilation.
%
%   `EnableImreadJpeg`:: `true`
%      Set this option to `true` to compile `vl_imreadjpeg`.
%
%   `EnableDouble`:: `true`
%      Set this option to `true` to compile the support for DOUBLE
%      data types.
%
%   `ImageLibrary`:: `libjpeg` (Linux), `gdiplus` (Windows), `quartz` (Mac)
%      The image library to use for `vl_impreadjpeg`.
%
%   `ImageLibraryCompileFlags`:: platform dependent
%      A cell-array of additional flags to use when compiling
%      `vl_imreadjpeg`.
%
%   `ImageLibraryLinkFlags`:: platform dependent
%      A cell-array of additional flags to use when linking
%      `vl_imreadjpeg`.
%
%   `EnableCudnn`:: `false`
%      Set to `true` to compile CuDNN support. See CuDNN
%      documentation for the Hardware/CUDA version requirements.
%
%   `CudnnRoot`:: `'local/'`
%      Directory containing the unpacked binaries and header files of
%      the CuDNN library.
%
%   ## Compiling the CPU code
%
%   By default, the `EnableGpu` option is switched to off, such that
%   the GPU code support is not compiled in.
%
%   Generally, you only need a 64bit C/C++ compiler (usually Xcode, GCC or
%   Visual Studio for Mac, Linux, and Windows respectively). The
%   compiler can be setup in MATLAB using the
%
%      mex -setup
%
%   command.
%
%   ## Compiling the GPU code
%
%   In order to compile the GPU code, set the `EnableGpu` option to
%   `true`. For this to work you will need:
%
%   * To satisfy all the requirements to compile the CPU code (see
%     above).
%
%   * A NVIDIA GPU with at least *compute capability 2.0*.
%
%   * The *MATALB Parallel Computing Toolbox*. This can be purchased
%     from Mathworks (type `ver` in MATLAB to see if this toolbox is
%     already comprised in your MATLAB installation; it often is).
%
%   * A copy of the *CUDA Devkit*, which can be downloaded for free
%     from NVIDIA. Note that each MATLAB version requires a
%     particular CUDA Devkit version:
%
%     | MATLAB version | Release | CUDA Devkit  |
%     |----------------|---------|--------------|
%     | 8.2            | 2013b   | 5.5          |
%     | 8.3            | 2014a   | 5.5          |
%     | 8.4            | 2014b   | 6.0          |
%     | 8.6            | 2015b   | 7.0          |
%     | 9.0            | 2016a   | 7.5          |
%
%     Different versions of CUDA may work using the hack described
%     above (i.e. setting the `CudaMethod` to `nvcc`).
%
%   The following configurations have been tested successfully:
%
%   * Windows 7 x64, MATLAB R2014a, Visual C++ 2010, 2013 and CUDA Toolkit
%     6.5. VS 2015 CPU version only (not supported by CUDA Toolkit yet).
%   * Windows 8 x64, MATLAB R2014a, Visual C++ 2013 and CUDA
%     Toolkit 6.5.
%   * Mac OS X 10.9, 10.10, 10.11, MATLAB R2013a to R2016a, Xcode, CUDA
%     Toolkit 5.5 to 7.5.
%   * GNU/Linux, MATALB R2014a/R2015a/R2015b/R2016a, gcc/g++, CUDA Toolkit 5.5/6.5/7.5.
%
%   Compilation on Windows with MinGW compiler (the default mex compiler in
%   Matlab) is not supported. For Windows, please reconfigure mex to use
%   Visual Studio C/C++ compiler.
%   Furthermore your GPU card must have ComputeCapability >= 2.0 (see
%   output of `gpuDevice()`) in order to be able to run the GPU code.
%   To change the compute capabilities, for `mex` `CudaMethod` edit
%   the particular config file.  For the 'nvcc' method, compute
%   capability is guessed based on the GPUDEVICE output. You can
%   override it by setting the 'CudaArch' parameter (e.g. in case of
%   multiple GPUs with various architectures).
%
%   See also: [Compliling MatConvNet](../install.md#compiling),
%   [Compiling MEX files containing CUDA
%   code](http://mathworks.com/help/distcomp/run-mex-functions-containing-cuda-code.html),
%   `vl_setup()`, `vl_imreadjpeg()`.

% Copyright (C) 2014-16 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% Get MatConvNet root directory
root = fileparts(fileparts(mfilename('fullpath'))) ;
addpath(fullfile(root, 'matlab')) ;

% --------------------------------------------------------------------
%                                                        Parse options
% --------------------------------------------------------------------

opts.enableGpu        = false;
opts.enableImreadJpeg = true;
opts.enableCudnn      = false;
opts.enableDouble     = true;
opts.imageLibrary = [] ;
opts.imageLibraryCompileFlags = {} ;
opts.imageLibraryLinkFlags = [] ;
opts.verbose          = 0;
opts.debug            = false;
opts.cudaMethod       = [] ;
opts.cudaRoot         = [] ;
opts.cudaArch         = [] ;
opts.defCudaArch      = [...
  '-gencode=arch=compute_20,code=\"sm_20,compute_20\" '...
  '-gencode=arch=compute_30,code=\"sm_30,compute_30\"'];
opts.cudnnRoot        = 'local/cudnn' ;
opts = vl_argparse(opts, varargin);

% --------------------------------------------------------------------
%                                                     Files to compile
% --------------------------------------------------------------------

arch = computer('arch') ;
if isempty(opts.imageLibrary)
  switch arch
    case 'glnxa64', opts.imageLibrary = 'libjpeg' ;
    case 'maci64', opts.imageLibrary = 'quartz' ;
    case 'win64', opts.imageLibrary = 'gdiplus' ;
  end
end
if isempty(opts.imageLibraryLinkFlags)
  switch opts.imageLibrary
    case 'libjpeg', opts.imageLibraryLinkFlags = {'-ljpeg'} ;
    case 'quartz', opts.imageLibraryLinkFlags = {'-framework Cocoa -framework ImageIO'} ;
    case 'gdiplus', opts.imageLibraryLinkFlags = {'gdiplus.lib'} ;
  end
end

lib_src = {} ;
mex_src = {} ;

% Files that are compiled as CPP or CU depending on whether GPU support
% is enabled.
if opts.enableGpu, ext = 'cu' ; else ext='cpp' ; end
lib_src{end+1} = fullfile(root,'matlab','src','bits',['data.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['datamex.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnconv.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnfullyconnected.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnsubsample.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnpooling.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnnormalize.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnbnorm.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnbias.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnbilinearsampler.' ext]) ;
lib_src{end+1} = fullfile(root,'matlab','src','bits',['nnroipooling.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnconv.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnconvt.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnpool.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnnormalize.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnbnorm.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnbilinearsampler.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_nnroipool.' ext]) ;
mex_src{end+1} = fullfile(root,'matlab','src',['vl_taccummex.' ext]) ;
switch arch
  case {'glnxa64','maci64'}
    % not yet supported in windows
    mex_src{end+1} = fullfile(root,'matlab','src',['vl_tmove.' ext]) ;
end

% CPU-specific files
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','im2row_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','subsample_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','copy_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','pooling_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','normalize_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','bnorm_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','tinythread.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','bilinearsampler_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','roipooling_cpu.cpp') ;
lib_src{end+1} = fullfile(root,'matlab','src','bits','imread.cpp') ;

% GPU-specific files
if opts.enableGpu
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','im2row_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','subsample_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','copy_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','pooling_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','normalize_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','bnorm_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','bilinearsampler_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','roipooling_gpu.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','datacu.cu') ;
  mex_src{end+1} = fullfile(root,'matlab','src','vl_cudatool.cu') ;
end

% cuDNN-specific files
if opts.enableCudnn
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','nnconv_cudnn.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','nnbias_cudnn.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','nnpooling_cudnn.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','nnbilinearsampler_cudnn.cu') ;
  lib_src{end+1} = fullfile(root,'matlab','src','bits','impl','nnbnorm_cudnn.cu') ;
end

% Other files
if opts.enableImreadJpeg
  mex_src{end+1} = fullfile(root,'matlab','src', ['vl_imreadjpeg.' ext]) ;
  mex_src{end+1} = fullfile(root,'matlab','src', ['vl_imreadjpeg_old.' ext]) ;
  lib_src{end+1} = fullfile(root,'matlab','src', 'bits', 'impl', ['imread_' opts.imageLibrary '.cpp']) ;
end

% --------------------------------------------------------------------
%                                                   Setup CUDA toolkit
% --------------------------------------------------------------------

if opts.enableGpu
  opts.verbose && fprintf('%s: * CUDA configuration *\n', mfilename) ;

  % Find the CUDA Devkit
  if isempty(opts.cudaRoot), opts.cudaRoot = search_cuda_devkit(opts) ; end
  opts.verbose && fprintf('%s:\tCUDA: using CUDA Devkit ''%s''.\n', ...
                          mfilename, opts.cudaRoot) ;

  opts.nvccPath = fullfile(opts.cudaRoot, 'bin', 'nvcc') ;
  switch arch
    case 'win64', opts.cudaLibDir = fullfile(opts.cudaRoot, 'lib', 'x64') ;
    case 'maci64', opts.cudaLibDir = fullfile(opts.cudaRoot, 'lib') ;
    case 'glnxa64', opts.cudaLibDir = fullfile(opts.cudaRoot, 'lib64') ;
    otherwise, error('Unsupported architecture ''%s''.', arch) ;
  end

  % Set the nvcc method as default for Win platforms
  if strcmp(arch, 'win64') && isempty(opts.cudaMethod)
    opts.cudaMethod = 'nvcc';
  end

  % Activate the CUDA Devkit
  cuver = activate_nvcc(opts.nvccPath) ;
  opts.verbose && fprintf('%s:\tCUDA: using NVCC ''%s'' (%d).\n', ...
                          mfilename, opts.nvccPath, cuver) ;

  % Set the CUDA arch string (select GPU architecture)
  if isempty(opts.cudaArch), opts.cudaArch = get_cuda_arch(opts) ; end
  opts.verbose && fprintf('%s:\tCUDA: NVCC architecture string: ''%s''.\n', ...
                          mfilename, opts.cudaArch) ;
end

if opts.enableCudnn
  opts.cudnnIncludeDir = fullfile(opts.cudnnRoot, 'include') ;
  switch arch
    case 'win64', opts.cudnnLibDir = fullfile(opts.cudnnRoot, 'lib', 'x64') ;
    case 'maci64', opts.cudnnLibDir = fullfile(opts.cudnnRoot, 'lib') ;
    case 'glnxa64', opts.cudnnLibDir = fullfile(opts.cudnnRoot, 'lib64') ;
    otherwise, error('Unsupported architecture ''%s''.', arch) ;
  end
end

% --------------------------------------------------------------------
%                                                     Compiler options
% --------------------------------------------------------------------

% Build directories
mex_dir = fullfile(root, 'matlab', 'mex') ;
bld_dir = fullfile(mex_dir, '.build');
if ~exist(fullfile(bld_dir,'bits','impl'), 'dir')
  mkdir(fullfile(bld_dir,'bits','impl')) ;
end

% Compiler flags
flags.cc = {} ;
flags.ccpass = {} ;
flags.ccoptim = {} ;
flags.link = {} ;
flags.linklibs = {} ;
flags.linkpass = {} ;
flags.nvccpass = {char(opts.cudaArch)} ;

if opts.verbose > 1
  flags.cc{end+1} = '-v' ;
end
if opts.debug
  flags.cc{end+1} = '-g' ;
  flags.nvccpass{end+1} = '-O0' ;
else
  flags.cc{end+1} = '-DNDEBUG' ;
  flags.nvccpass{end+1} = '-O3' ;
end
if opts.enableGpu
  flags.cc{end+1} = '-DENABLE_GPU' ;
end
if opts.enableCudnn
  flags.cc{end+1} = '-DENABLE_CUDNN' ;
  flags.cc{end+1} = ['-I"' opts.cudnnIncludeDir '"'] ;
end
if opts.enableDouble
  flags.cc{end+1} = '-DENABLE_DOUBLE' ;
end
flags.link{end+1} = '-lmwblas' ;
switch arch
  case {'maci64'}
  case {'glnxa64'}
    flags.linklibs{end+1} = '-lrt' ;
  case {'win64'}
    % VisualC does not pass this even if available in the CPU architecture
    flags.cc{end+1} = '-D__SSSE3__' ;
end

if opts.enableImreadJpeg
  flags.cc = horzcat(flags.cc, opts.imageLibraryCompileFlags) ;
  flags.linklibs = horzcat(flags.linklibs, opts.imageLibraryLinkFlags) ;
end

if opts.enableGpu
  flags.link = horzcat(flags.link, {['-L"' opts.cudaLibDir '"'], '-lcudart', '-lcublas'}) ;
  switch arch
    case {'maci64', 'glnxa64'}
      flags.link{end+1} = '-lmwgpu' ;
    case 'win64'
      flags.link{end+1} = '-lgpu' ;
  end
  if opts.enableCudnn
    flags.link{end+1} = ['-L"' opts.cudnnLibDir '"'] ;
    flags.link{end+1} = '-lcudnn' ;
  end
end

switch arch
  case {'maci64'}
    flags.ccpass{end+1} = '-mmacosx-version-min=10.9' ;
    flags.linkpass{end+1} = '-mmacosx-version-min=10.9' ;
    flags.ccoptim{end+1} = '-mssse3 -ffast-math' ;
    flags.nvccpass{end+1} = '-Xcompiler -fPIC' ;

    if opts.enableGpu
      flags.linkpass{end+1} = sprintf('-Wl,-rpath -Wl,"%s"', opts.cudaLibDir) ;
    end
    if opts.enableGpu && opts.enableCudnn
      flags.linkpass{end+1} = sprintf('-Wl,-rpath -Wl,"%s"', opts.cudnnLibDir) ;
    end

    if opts.enableGpu && cuver < 70000
      % CUDA prior to 7.0 on Mac require GCC libstdc++ instead of the native
      % clang libc++. This should go away in the future.
      flags.ccpass{end+1} = '-stdlib=libstdc++' ;
      flags.linkpass{end+1} = '-stdlib=libstdc++' ;
      if  ~verLessThan('matlab', '8.5.0')
        % Complicating matters, MATLAB 8.5.0 links to clang's libc++ by
        % default when linking MEX files overriding the option above. We
        % force it to use GCC libstdc++
        flags.linkpass{end+1} = '-L"$MATLABROOT/bin/maci64" -lmx -lmex -lmat -lstdc++' ;
      end
    end

  case {'glnxa64'}
    flags.ccoptim{end+1} = '-mssse3 -ftree-vect-loop-version -ffast-math -funroll-all-loops' ;
    flags.nvccpass{end+1} = '-Xcompiler -fPIC -D_FORCE_INLINES' ;

    if opts.enableGpu
      flags.linkpass{end+1} = sprintf('-Wl,-rpath -Wl,"%s"', opts.cudaLibDir) ;
    end
    if opts.enableGpu && opts.enableCudnn
      flags.linkpass{end+1} = sprintf('-Wl,-rpath -Wl,"%s"', opts.cudnnLibDir) ;
    end

  case {'win64'}
    flags.nvccpass{end+1} = '-Xcompiler /MD' ;
    cl_path = fileparts(check_clpath()); % check whether cl.exe in path
    flags.nvccpass{end+1} = sprintf('--compiler-bindir "%s"', cl_path) ;
end

% --------------------------------------------------------------------
%                                                        Command flags
% --------------------------------------------------------------------

flags.mexcc = horzcat(flags.cc, ...
                      {'-largeArrayDims'}, ...
                      {['CXXFLAGS=$CXXFLAGS ' strjoin(flags.ccpass)]}, ...
                      {['CXXOPTIMFLAGS=$CXXOPTIMFLAGS ' strjoin(flags.ccoptim)]}) ;
if ~ispc, flags.mexcc{end+1} = '-cxx'; end

% mex: compile GPU
flags.mexcu= horzcat({'-f' mex_cuda_config(root)}, ...
                     flags.cc, ...
                     {'-largeArrayDims'}, ...
                     {['CXXFLAGS=$CXXFLAGS ' quote_nvcc(flags.ccpass) ' ' strjoin(flags.nvccpass)]}, ...
                     {['CXXOPTIMFLAGS=$CXXOPTIMFLAGS ' quote_nvcc(flags.ccoptim)]}) ;

% mex: link
flags.mexlink = horzcat(flags.cc, flags.link, ...
                        {'-largeArrayDims'}, ...
                        {['LDFLAGS=$LDFLAGS ', strjoin(flags.linkpass)]}, ...
                        {['LINKLIBS=', strjoin(flags.linklibs), ' $LINKLIBS']}) ;

% nvcc: compile GPU
flags.nvcc = horzcat(flags.cc, ...
                     {opts.cudaArch}, ...
                     {sprintf('-I"%s"', fullfile(matlabroot, 'extern', 'include'))}, ...
                     {sprintf('-I"%s"', fullfile(matlabroot, 'toolbox','distcomp','gpu','extern','include'))}, ...
                     {quote_nvcc(flags.ccpass)}, ...
                     {quote_nvcc(flags.ccoptim)}, ...
                     flags.nvccpass) ;

if opts.verbose
  fprintf('%s: * Compiler and linker configurations *\n', mfilename) ;
  fprintf('%s: \tintermediate build products directory: %s\n', mfilename, bld_dir) ;
  fprintf('%s: \tMEX files: %s/\n', mfilename, mex_dir) ;
  fprintf('%s: \tMEX options [CC CPU]: %s\n', mfilename, strjoin(flags.mexcc)) ;
  fprintf('%s: \tMEX options [LINK]: %s\n', mfilename, strjoin(flags.mexlink)) ;
end
if opts.verbose && opts.enableGpu
  fprintf('%s: \tMEX options [CC GPU]: %s\n', mfilename, strjoin(flags.mexcu)) ;
end
if opts.verbose && opts.enableGpu && strcmp(opts.cudaMethod,'nvcc')
  fprintf('%s: \tNVCC options [CC GPU]: %s\n', mfilename, strjoin(flags.nvcc)) ;
end
if opts.verbose && opts.enableImreadJpeg
  fprintf('%s: * Reading images *\n', mfilename) ;
  fprintf('%s: \tvl_imreadjpeg enabled\n', mfilename) ;
  fprintf('%s: \timage library: %s\n', mfilename, opts.imageLibrary) ;
  fprintf('%s: \timage library compile flags: %s\n', mfilename, strjoin(opts.imageLibraryCompileFlags)) ;
  fprintf('%s: \timage library link flags: %s\n', mfilename, strjoin(opts.imageLibraryLinkFlags)) ;
end

% --------------------------------------------------------------------
%                                                              Compile
% --------------------------------------------------------------------

% Intermediate object files
srcs = horzcat(lib_src,mex_src) ;
for i = 1:numel(horzcat(lib_src, mex_src))
  [~,~,ext] = fileparts(srcs{i}) ; ext(1) = [] ;
  objfile = toobj(bld_dir,srcs{i});
  if strcmp(ext,'cu')
    if strcmp(opts.cudaMethod,'nvcc')
      nvcc_compile(opts, srcs{i}, objfile, flags.nvcc) ;
    else
      mex_compile(opts, srcs{i}, objfile, flags.mexcu) ;
    end
  else
    mex_compile(opts, srcs{i}, objfile, flags.mexcc) ;
  end
  assert(exist(objfile, 'file') ~= 0, 'Compilation of %s failed.', objfile);
end

% Link into MEX files
for i = 1:numel(mex_src)
  objs = toobj(bld_dir, [mex_src(i), lib_src]) ;
  mex_link(opts, objs, mex_dir, flags.mexlink) ;
end

% Reset path adding the mex subdirectory just created
vl_setupnn() ;

% --------------------------------------------------------------------
%                                                    Utility functions
% --------------------------------------------------------------------

% --------------------------------------------------------------------
function objs = toobj(bld_dir,srcs)
% --------------------------------------------------------------------
str = fullfile('matlab','src') ;
multiple = iscell(srcs) ;
if ~multiple, srcs = {srcs} ; end
objs = cell(1, numel(srcs));
for t = 1:numel(srcs)
  i = strfind(srcs{t},str);
  objs{t} = fullfile(bld_dir, srcs{t}(i+numel(str):end)) ;
end
if ~multiple, objs = objs{1} ; end
objs = regexprep(objs,'.cpp$',['.' objext]) ;
objs = regexprep(objs,'.cu$',['.' objext]) ;
objs = regexprep(objs,'.c$',['.' objext]) ;

% --------------------------------------------------------------------
function mex_compile(opts, src, tgt, mex_opts)
% --------------------------------------------------------------------
mopts = {'-outdir', fileparts(tgt), src, '-c', mex_opts{:}} ;
opts.verbose && fprintf('%s: MEX CC: %s\n', mfilename, strjoin(mopts)) ;
mex(mopts{:}) ;

% --------------------------------------------------------------------
function nvcc_compile(opts, src, tgt, nvcc_opts)
% --------------------------------------------------------------------
nvcc_path = fullfile(opts.cudaRoot, 'bin', 'nvcc');
nvcc_cmd = sprintf('"%s" -c "%s" %s -o "%s"', ...
                   nvcc_path, src, ...
                   strjoin(nvcc_opts), tgt);
opts.verbose && fprintf('%s: NVCC CC: %s\n', mfilename, nvcc_cmd) ;
status = system(nvcc_cmd);
if status, error('Command %s failed.', nvcc_cmd); end;

% --------------------------------------------------------------------
function mex_link(opts, objs, mex_dir, mex_flags)
% --------------------------------------------------------------------
mopts = {'-outdir', mex_dir, mex_flags{:}, objs{:}} ;
opts.verbose && fprintf('%s: MEX LINK: %s\n', mfilename, strjoin(mopts)) ;
mex(mopts{:}) ;

% --------------------------------------------------------------------
function ext = objext()
% --------------------------------------------------------------------
% Get the extension for an 'object' file for the current computer
% architecture
switch computer('arch')
  case 'win64', ext = 'obj';
  case {'maci64', 'glnxa64'}, ext = 'o' ;
  otherwise, error('Unsupported architecture %s.', computer) ;
end

% --------------------------------------------------------------------
function conf_file = mex_cuda_config(root)
% --------------------------------------------------------------------
% Get mex CUDA config file
mver = [1e4 1e2 1] * sscanf(version, '%d.%d.%d') ;
if mver <= 80200, ext = 'sh' ; else ext = 'xml' ; end
arch = computer('arch') ;
switch arch
  case {'win64'}
    config_dir = fullfile(matlabroot, 'toolbox', ...
                          'distcomp', 'gpu', 'extern', ...
                          'src', 'mex', arch) ;
  case {'maci64', 'glnxa64'}
    config_dir = fullfile(root, 'matlab', 'src', 'config') ;
end
conf_file = fullfile(config_dir, ['mex_CUDA_' arch '.' ext]);
fprintf('%s:\tCUDA: MEX config file: ''%s''\n', mfilename, conf_file);

% --------------------------------------------------------------------
function cl_path = check_clpath()
% --------------------------------------------------------------------
% Checks whether the cl.exe is in the path (needed for the nvcc). If
% not, tries to guess the location out of mex configuration.
cc = mex.getCompilerConfigurations('c++');
if isempty(cc)
  error(['Mex is not configured.'...
    'Run "mex -setup" to configure your compiler. See ',...
    'http://www.mathworks.com/support/compilers ', ...
    'for supported compilers for your platform.']);
end
cl_path = fullfile(cc.Location, 'VC', 'bin', 'amd64');
[status, ~] = system('cl.exe -help');
if status == 1
  warning('CL.EXE not found in PATH. Trying to guess out of mex setup.');
  prev_path = getenv('PATH');
  setenv('PATH', [prev_path ';' cl_path]);
  status = system('cl.exe');
  if status == 1
    setenv('PATH', prev_path);
    error('Unable to find cl.exe');
  else
    fprintf('Location of cl.exe (%s) successfully added to your PATH.\n', ...
      cl_path);
  end
end

% -------------------------------------------------------------------------
function paths = which_nvcc()
% -------------------------------------------------------------------------
switch computer('arch')
  case 'win64'
    [~, paths] = system('where nvcc.exe');
    paths = strtrim(paths);
    paths = paths(strfind(paths, '.exe'));
  case {'maci64', 'glnxa64'}
    [~, paths] = system('which nvcc');
    paths = strtrim(paths) ;
end

% -------------------------------------------------------------------------
function cuda_root = search_cuda_devkit(opts)
% -------------------------------------------------------------------------
% This function tries to to locate a working copy of the CUDA Devkit.

opts.verbose && fprintf(['%s:\tCUDA: searching for the CUDA Devkit' ...
                    ' (use the option ''CudaRoot'' to override):\n'], mfilename);

% Propose a number of candidate paths for NVCC
paths = {getenv('MW_NVCC_PATH')} ;
paths = [paths, which_nvcc()] ;
for v = {'5.5', '6.0', '6.5', '7.0', '7.5', '8.0', '8.5', '9.0'}
  switch computer('arch')
    case 'glnxa64'
      paths{end+1} = sprintf('/usr/local/cuda-%s/bin/nvcc', char(v)) ;
    case 'maci64'
      paths{end+1} = sprintf('/Developer/NVIDIA/CUDA-%s/bin/nvcc', char(v)) ;
    case 'win64'
      paths{end+1} = sprintf('C:\\Program Files\\NVIDIA GPU Computing Toolkit\\CUDA\\v%s\\bin\\nvcc.exe', char(v)) ;
  end
end
paths{end+1} = sprintf('/usr/local/cuda/bin/nvcc') ;

% Validate each candidate NVCC path
for i=1:numel(paths)
  nvcc(i).path = paths{i} ;
  [nvcc(i).isvalid, nvcc(i).version] = validate_nvcc(paths{i}) ;
end
if opts.verbose
  fprintf('\t| %5s | %5s | %-70s |\n', 'valid', 'ver', 'NVCC path') ;
  for i=1:numel(paths)
    fprintf('\t| %5d | %5d | %-70s |\n', ...
            nvcc(i).isvalid, nvcc(i).version, nvcc(i).path) ;
  end
end

% Pick an entry
index = find([nvcc.isvalid]) ;
if isempty(index)
  error('Could not find a valid NVCC executable\n') ;
end
[~, newest] = max([nvcc(index).version]);
nvcc = nvcc(index(newest)) ;
cuda_root = fileparts(fileparts(nvcc.path)) ;

if opts.verbose
  fprintf('%s:\tCUDA: choosing NVCC compiler ''%s'' (version %d)\n', ...
          mfilename, nvcc.path, nvcc.version) ;
end

% -------------------------------------------------------------------------
function [valid, cuver]  = validate_nvcc(nvccPath)
% -------------------------------------------------------------------------
[status, output] = system(sprintf('"%s" --version', nvccPath)) ;
valid = (status == 0) ;
if ~valid
  cuver = 0 ;
  return ;
end
match = regexp(output, 'V(\d+\.\d+\.\d+)', 'match') ;
if isempty(match), valid = false ; return ; end
cuver = [1e4 1e2 1] * sscanf(match{1}, 'V%d.%d.%d') ;

% --------------------------------------------------------------------
function cuver = activate_nvcc(nvccPath)
% --------------------------------------------------------------------

% Validate the NVCC compiler installation
[valid, cuver] = validate_nvcc(nvccPath) ;
if ~valid
  error('The NVCC compiler ''%s'' does not appear to be valid.', nvccPath) ;
end

% Make sure that NVCC is visible by MEX by setting the MW_NVCC_PATH
% environment variable to the NVCC compiler path
if ~strcmp(getenv('MW_NVCC_PATH'), nvccPath)
  warning('Setting the ''MW_NVCC_PATH'' environment variable to ''%s''', nvccPath) ;
  setenv('MW_NVCC_PATH', nvccPath) ;
end

% In some operating systems and MATLAB versions, NVCC must also be
% available in the command line search path. Make sure that this is%
% the case.
[valid_, cuver_] = validate_nvcc('nvcc') ;
if ~valid_ || cuver_ ~= cuver
  warning('NVCC not found in the command line path or the one found does not matches ''%s''.', nvccPath);
  nvccDir = fileparts(nvccPath) ;
  prevPath = getenv('PATH') ;
  switch computer
    case 'PCWIN64', separator = ';' ;
    case {'GLNXA64', 'MACI64'}, separator = ':' ;
  end
  setenv('PATH', [nvccDir separator prevPath]) ;
  [valid_, cuver_] = validate_nvcc('nvcc') ;
  if ~valid_ || cuver_ ~= cuver
    setenv('PATH', prevPath) ;
    error('Unable to set the command line path to point to ''%s'' correctly.', nvccPath) ;
  else
    fprintf('Location of NVCC (%s) added to your command search PATH.\n', nvccDir) ;
  end
end

% -------------------------------------------------------------------------
function str = quote_nvcc(str)
% -------------------------------------------------------------------------
if iscell(str), str = strjoin(str) ; end
str = strrep(strtrim(str), ' ', ',') ;
if ~isempty(str), str = ['-Xcompiler ' str] ; end

% --------------------------------------------------------------------
function cudaArch = get_cuda_arch(opts)
% --------------------------------------------------------------------
opts.verbose && fprintf('%s:\tCUDA: determining GPU compute capability (use the ''CudaArch'' option to override)\n', mfilename);
try
  gpu_device = gpuDevice();
  arch_code = strrep(gpu_device.ComputeCapability, '.', '');
  cudaArch = ...
      sprintf('-gencode=arch=compute_%s,code=\\\"sm_%s,compute_%s\\\" ', ...
              arch_code, arch_code, arch_code) ;
catch
  opts.verbose && fprintf(['%s:\tCUDA: cannot determine the capabilities of the installed GPU; ' ...
                      'falling back to default\n'], mfilename);
  cudaArch = opts.defCudaArch;
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function camPaths = dtiLoadCaminoPaths(filename, xform, bEndOnly)
%
% camPaths = dtiLoadCaminoPaths(filename)
%
% Loads fiber pathways from a Camino .Bfloat file.
%
% Notes: Pathways are assumed to be stored in the following Camino format:
% [<N>, <seed point index>, <x_1>, <y_1>, <z_1>,...,<x_numPoints>, <y_N>,
% <z_N>, <N>,...,<z_N>]
%
%
% HISTORY:
% 2007.09.07 Written by Anthony Sherbondy
% 2010.04.19 AJS made the code much faster by removing function calls and
% preallocating cell array.

% %% Endpoints only?
% if ieNotDefined('bEndOnly')
%     bEndOnly = 0;
% end

%% Open file and read all raw data in
fid = fopen(filename,'rb','b');
rawPaths = fread(fid,'float'); fclose(fid);

%% Create pathway structure
camPaths.statHeader(1).agg_name = 'Length';
camPaths.pathwayInfo = [];
camPaths.pathways = {};

%% First read off all pathway sizes
buffer_size = length(rawPaths);
v_num_points = nan(buffer_size,1);
v_sptr = zeros(buffer_size,1);
sPtr = 1;
pp = 1;
while sPtr < buffer_size
    v_num_points(pp) = rawPaths(sPtr);
    v_sptr(pp) = (sPtr);
    sPtr = sPtr+3*v_num_points(pp)+2;
    pp = pp + 1;
end

%% Get acutal number of pathways represented in buffers
v_sptr = v_sptr(~isnan(v_num_points) & v_num_points>=3);
v_num_points = v_num_points(~isnan(v_num_points) & v_num_points>=3);

%% Now lets go and read all the pathways
camPaths.pathways = cell(numel(v_num_points),1);
%pct_step = 10;
%pct = pct_step;
%step = floor(length(v_num_points)/(100/pct_step));
for pp=1:length(v_num_points)
    sPtr = v_sptr(pp);
    numPoints = v_num_points(pp);
    camPaths.pathways{pp,1} = reshape(rawPaths(sPtr+2:sPtr+1+3*numPoints),3,numPoints);
%     if pp > step * pct
%         pct = floor(pp / step);
%         fprintf(1,'\n%3d percent file processed...', pct);
%         pct = pct + pct_step;
%     end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = process_spikesorting_ultramegasort2000( varargin )
% PROCESS_SPIKESORTING_ULTRAMEGASORT2000:
% This process separates the initial raw signal to nChannels binary signals
% and performs spike sorting individually on each channel with the
% UltraMegaSort2000 spike-sorter. The spikes are clustered and assigned to
% individual neurons. The code ultimately produces a raw_elec(i)_spikes.mat
% for each electrode that can be used later for supervised spike-sorting.
% When all spikes on all electrodes have been clustered, all the spikes for
% each neuron is assigned to an events file in brainstorm format.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Konstantinos Nasiotis, 2018-2022
%          Martin Cousineau, 2018
%          Francois Tadel, 2022

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'UltraMegaSort2000';
    sProcess.Category    = 'File';
    sProcess.SubGroup    = {'Electrophysiology','Unsupervised Spike Sorting'};
    sProcess.Index       = 1202;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/e-phys/SpikeSorting';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw'};
    sProcess.OutputTypes = {'raw'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 0;
    % Spike sorter name
    sProcess.options.spikesorter.Type   = 'text';
    sProcess.options.spikesorter.Value  = 'ultramegasort2000';
    sProcess.options.spikesorter.Hidden = 1;
    % RAM limitation
    sProcess.options.binsize.Comment = 'Maximum RAM to use: ';
    sProcess.options.binsize.Type    = 'value';
    sProcess.options.binsize.Value   = {2, 'GB', 1};
    % Parallel processing
    sProcess.options.parallel.Comment = 'Parallel processing';
    sProcess.options.parallel.Type    = 'checkbox';
    sProcess.options.parallel.Value   = 0;
    % Use SSP/ICA
    sProcess.options.usessp.Comment = 'Apply the existing SSP/ICA projectors';
    sProcess.options.usessp.Type    = 'checkbox';
    sProcess.options.usessp.Value   = 1;
    % Separator
    sProcess.options.sep1.Type = 'label';
    sProcess.options.sep1.Comment = '
';
    % === Low bound
    sProcess.options.highpass.Comment = 'Lower cutoff frequency:';
    sProcess.options.highpass.Type    = 'value';
    sProcess.options.highpass.Value   = {700,'Hz ',0};
    % === High bound
    sProcess.options.lowpass.Comment = 'Upper cutoff frequency:';
    sProcess.options.lowpass.Type    = 'value';
    sProcess.options.lowpass.Value   = {5000,'Hz ',0};
    % Separator
    sProcess.options.sep2.Type = 'label';
    sProcess.options.sep2.Comment = '
';
    % Options: Options
    sProcess.options.edit.Comment = {'panel_spikesorting_options', 'UltraMegaSort2000 parameters: '};
    sProcess.options.edit.Type    = 'editpref';
    sProcess.options.edit.Value   = [];
    % Label: Reset options
    sProcess.options.edit_help.Comment = 'To restore default options: re-install the ultramegasort plugin.';
    sProcess.options.edit_help.Type    = 'label';
    % Label: Warning that pre-spikesorted events will be overwritten
    sProcess.options.warning.Comment = '
Warning: Existing spike events will be overwritten';
    sProcess.options.warning.Type    = 'label';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInput) %#ok<DEFNU>
    OutputFiles = {};

    % ===== DEPENDENCIES =====
    % Not available in the compiled version
    if bst_iscompiled()
        bst_report('Error', sProcess, sInput, 'This function is not available in the compiled version of Brainstorm.');
        return
    end
    % Load plugin
    [isInstalled, errMsg] = bst_plugin('Install', 'ultramegasort2000');
    if ~isInstalled
        error(errMsg);
    end

    % ===== OPTIONS =====
    % Get option: bin size
    BinSize = sProcess.options.binsize.Value{1};
    if (BinSize <= 0)
        bst_report('Error', sProcess, sInput, 'Invalid maximum amount of RAM specified.');
        return
    end
    % Get other options
    isParallel = sProcess.options.parallel.Value;
    UseSsp = sProcess.options.usessp.Value;
    LowPass = sProcess.options.lowpass.Value{1}(1);
    HighPass = sProcess.options.highpass.Value{1}(1);

    % ===== LOAD INPUTS =====
    % Get protocol info
    ProtocolInfo = bst_get('ProtocolInfo');
    BrainstormTmpDir = bst_get('BrainstormTmpDir');
    % File path
    [fPath, fBase] = bst_fileparts(file_fullpath(sInput.FileName));
    % Remove "data_0raw" or "data_" tag
    if (length(fBase) > 10 && strcmp(fBase(1:10), 'data_0raw_'))
        fBase = fBase(11:end);
    elseif (length(fBase) > 5) && strcmp(fBase(1:5), 'data_')
        fBase = fBase(6:end);
    end

    % Load input files
    DataMat = in_bst_data(sInput.FileName, 'F');
    sFile = DataMat.F;
    % Check filtering frequencies
    nyq = floor(sFile.prop.sfreq/2);
    if (LowPass >= nyq)
        bst_report('Error', sProcess, sInput, ['Higher cutoff frequency must be lower than Nyquist frequency (' num2str(nyq) ' Hz).']);
        return;
    elseif (HighPass >= LowPass)
        bst_report('Error', sProcess, sInput, 'Higher cutoff frequency must be lower lower cutoff frequency.');
        return;
    end
    % Load channel file
    ChannelMat = in_bst_channel(sInput.ChannelFile);
    numChannels = length(ChannelMat.Channel);
    % Demultiplex channels
    demultiplexDir = bst_fullfile(BrainstormTmpDir, 'Unsupervised_Spike_Sorting', ProtocolInfo.Comment, sInput.FileName);
    elecFiles = out_demultiplex(sInput.FileName, sInput.ChannelFile, demultiplexDir, UseSsp, BinSize * 1e9, isParallel);

    % ===== OUTPUT FOLDER =====
    outputPath = bst_fullfile(fPath, [fBase '_ums2k_spikes']);
    previous_directory = pwd;
    % If output folder already exists: delete it
    if exist(outputPath, 'dir') == 7
        % Move Matlab out of the folder to be deleted
        if ~isempty(strfind(previous_directory, outputPath))
            cd(bst_fileparts(outputPath));
        end
        % Delete existing output folder
        try
            rmdir(outputPath, 's');
        catch
        	error(['Could not remove spikes folder: ' 10 outputPath 10 ' Make sure this folder is not open in another program.'])
        end
    end
    % Create output folder
    mkdir(outputPath);

    % ===== SPIKE SORTING =====
    cd(outputPath);
    if isParallel
        bst_progress('start', 'Spike-sorting', 'Extracting spikes...');
        parfor ielectrode = 1:numChannels
            do_UltraMegaSorting(elecFiles{ielectrode}, sFile, LowPass, HighPass, sFile.prop.sfreq);
        end
    else
        bst_progress('start', 'Spike-sorting', 'Extracting spikes...', 0, numChannels);
        for ielectrode = 1:numChannels
            do_UltraMegaSorting(elecFiles{ielectrode}, sFile, LowPass, HighPass, sFile.prop.sfreq);
            bst_progress('inc', 1);
        end
    end
    % Restore current folder
    cd(previous_directory);

    % ===== IMPORT EVENTS =====
    bst_progress('start', 'UltraMegaSort2000', 'Gathering spiking events...');
    % Delete existing spike events
    panel_spikes('DeleteSpikeEvents', sInput.FileName);

    % Build output filename
    NewBstFilePrefix = bst_fullfile(fPath, ['data_0ephys_ums2k_' fBase]);
    NewBstFile = [NewBstFilePrefix '.mat'];
    iFile = 1;
    commentSuffix = '';
    while exist(NewBstFile, 'file') == 2
        iFile = iFile + 1;
        NewBstFile = [NewBstFilePrefix '_' num2str(iFile) '.mat'];
        commentSuffix = [' (' num2str(iFile) ')'];
    end
    % Build output structure
    DataMat = struct();
    DataMat.Comment     = ['UltraMegaSort2000 Spike Sorting' commentSuffix];
    DataMat.DataType    = 'raw';
    DataMat.Device      = 'ultramegasort2000';
    DataMat.Name        = file_short(NewBstFile);
    DataMat.Parent      = file_short(outputPath);
    DataMat.RawFile     = sInput.FileName;
    DataMat.Spikes      = struct();
    % Build spikes structure
    spikes = dir(bst_fullfile(outputPath, 'times_raw_elec*.mat'));
    spikes = sort_nat({spikes.name});
    for iSpike = 1:length(spikes)
        DataMat.Spikes(iSpike).Path = file_short(outputPath);
        DataMat.Spikes(iSpike).File = spikes{iSpike};
        if exist(bst_fullfile(outputPath, DataMat.Spikes(iSpike).File), 'file') ~= 2
            DataMat.Spikes(iSpike).File = '';
        end
        DataMat.Spikes(iSpike).Name = ChannelMat.Channel(iSpike).Name;
        DataMat.Spikes(iSpike).Mod  = 0;
    end

    % Save events file for backup
    SaveBrainstormEvents(DataMat, 'events_UNSUPERVISED.mat');

    % Add history field
    DataMat = bst_history('add', DataMat, 'import', ['Link to unsupervised electrophysiology files: ' outputPath]);
    % Save file on hard drive
    bst_save(NewBstFile, DataMat, 'v6');
    % Add file to database
    db_add_data(sInput.iStudy, file_short(NewBstFile), DataMat);
    % Return new file
    OutputFiles{end+1} = NewBstFile;

    % ===== UPDATE DATABASE =====
    % Update links
    db_links('Study', sInput.iStudy);
    panel_protocols('UpdateNode', 'Study', sInput.iStudy);
end

%% ===== SAVE BRAINSTORM EVENTS =====
function SaveBrainstormEvents(SpikeMat, outputFile, eventNamePrefix)
    if nargin < 3
        eventNamePrefix = '';
    end

    numElectrodes = length(SpikeMat.Spikes);
    iNewEvent = 0;
    events = struct();

    % Add existing non-spike events for backup
    DataMat = in_bst_data(SpikeMat.RawFile);
    existingEvents = DataMat.F.events;
    for iEvent = 1:length(existingEvents)
        if ~panel_spikes('IsSpikeEvent', existingEvents(iEvent).label)
            if iNewEvent == 0
                events = existingEvents(iEvent);
            else
                events(iNewEvent + 1) = existingEvents(iEvent);
            end
            iNewEvent = iNewEvent + 1;
        end
    end

    for iElectrode = 1:numElectrodes
        newEvents = panel_spikes(...
            'CreateSpikeEvents', ...
            SpikeMat.RawFile, ...
            SpikeMat.Device, ...
            bst_fullfile(file_fullpath(SpikeMat.Parent), SpikeMat.Spikes(iElectrode).File), ...
            SpikeMat.Spikes(iElectrode).Name, ...
            1, eventNamePrefix);

        if iNewEvent == 0
            events = newEvents;
            iNewEvent = length(newEvents);
        else
            numNewEvents = length(newEvents);
            events(iNewEvent+1:iNewEvent+numNewEvents) = newEvents;
            iNewEvent = iNewEvent + numNewEvents;
        end
    end

    save(bst_fullfile(file_fullpath(SpikeMat.Parent), outputFile),'events');
end

%% ===== ULTRAMEGA SORTING =====
function do_UltraMegaSorting(electrodeFile, sFile, LowPass, HighPass, Fs)
    [path, filename] = fileparts(electrodeFile);
    % Apply BST bandpass filter
    try
        DataMat = load(electrodeFile, 'data');
        filtered_data = bst_bandpass_hfilter(DataMat.data', Fs, HighPass, LowPass, 0, 0);
    catch e
        error(['Frequency filtering failed, try with a different frequency band.' 10 'Error: ' e.message]);
    end
    % Convert to a column vector in a cell array
    filtered_data = {filtered_data'};
    % Run spike sorting
    try
        spikes = ss_default_params(sFile.prop.sfreq);
        spikes = ss_detect(filtered_data,spikes);
        spikes = ss_align(spikes);
        spikes = ss_kmeans(spikes);
        spikes = ss_energy(spikes);
        spikes = ss_aggregate(spikes);
        save(['times_' filename '.mat'], 'spikes')
    catch e
        % If an error occurs, just don't create the spike file.
        clean_label = strrep(filename,'raw_elec_','');
        disp(e);
        disp(['Warning: Spiking failed on electrode ' clean_label '. Skipping this electrode.']);
    end
end

author: brainstorm-tools


text:
function fgOut=dtiGetFibersConnectingGM(fg, dt, minDist, showFig)
%Select fibers with both endpoints located near grey matter.
%
%  fg=dtiFibersConnectingGM(fg, dt, [minDist=1.74])
%
% Returns a fibers structure which is a subset of fibers from the original
% fiber group, keeping only the fibers whose both endpoints terminate near
% grey matter or the most inferior axial plane (the latter to trap
% cortico-spinal fibers). Relies on spm to perform tissue segmentation of b0.
% Grey matter mask is rather generous: anything with probability(gm)>0.
%
% Input parameters:
% minDist  - maximum distance from a fiber endpoint to Gm mask.
% showFig  - 'true', will show Gm mask; 'false', will save the mask.
%
% HISTORY:
% 07/31/2009 ER wrote it

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

if ~exist('minDist','var') || isempty(minDist)
    minDist=1.74;
end

[wm, gm, csf] = mrAnatSpmSegment(dt.b0, dt.xformToAcpc, 'mniepi'); gm=gm>=127;

% Only fibers whose both endpoints are within the gray matter (roi prepared) will be considered.
[x1,y1,z1] = ind2sub(size(gm), find(gm));

%fill up the most inferior nonzero slice & below with "gray matter" voxels -- to make
%sure corticospinal fibers will not be eliminated when retaining only
%"gray matter connecting" voxels.
gm_withcst=gm;
gm_withcst(:, :, min(z1))=1;

[x1_withcst,y1_withcst,z1_withcst] = ind2sub(size(gm_withcst), find(gm_withcst));
roi4cst= dtiNewRoi('mrAnatSpmSegment_gm');
roi4cst.coords = mrAnatXformCoords(dt.xformToAcpc, [x1_withcst,y1_withcst,z1_withcst]);

[fgOut] = dtiIntersectFibersWithRoi([], {'and', 'both_endpoints'}, minDist, roi4cst, fg);

if showFig
    showMontage(gm_withcst);
else
    imwrite(makeMontage(gm_withcst), fullfile(fileparts(dt.dataFile), 'bin', 'gmMask.png'));
end


text:
function [lats, lons] = GetCoastalPoints(Quality,maxLevel,minArea,LatLon)
%GETCOASTALPOINTS Reads the coastal points from a GSHHS file (version 2.x)
%
% USAGE:
%   [lats, lons] = GetCoastalPoints(Quality,maxLevel,minArea,LatLon)
%
% INPUTS:
%   Quality - integer specifying the quality of the coastline to load.
%       Acceptable values are 1-5. Only the coastline corresponding to 1
%       is available by default. REQUIRED
%   maxLevel - max level of data coastline to load. Levels are:
%           1 land, 2 lake, 3 island in lake, 4 pond in island in lake.
%           Default 1
%   minArea - minimum area (km^2) within a coastline to load. Default 10
%   LatLon - 1 by 4 vector describing a window in which a coastline must be
%       included in order for it to be loaded. Values are
%       [minLat minLong maxLat maxLon]. If empty load all. Default []
%
% OUTPUTS:
%   lats - latitude coordinates.
%   lons - longitude coordinates.
%
% Note: shapes are separated by NaN
%
% GSHHS data and a description of its format can be found at:
%
% "A Global Self-consistent, Hierarchical, High-resolution Shoreline Database"
% http://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html
%
% or
%
% "A Global Self-consistent, Hierarchical, High-resolution Geography Database"
% http://www.soest.hawaii.edu/pwessel/gshhs/index.html
%
% Author: Sean Hatch
% Modified by Wade Schwartzkopf to handle GSHHS version 2.x
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%Modified by Pat Cutler (NGA) 20130306 for load images segments LatLon
%input and documented inputs and outputs

if ~exist('maxLevel','var') || isempty(maxLevel)
    maxLevel = 1;
end

if ~exist('minArea','var') || isempty(minArea)
    minArea = 10;
end

if ~exist('LatLon','var') || isempty(LatLon)
    LatLonCheck = 0;
else
    LatLonCheck = 1;
end

% open the GSHHS data file
% the _c file contains the crudest data in the shoreline hierarchy
files = {'gshhs_c.b', 'gshhs_l.b', 'gshhs_i.b', 'gshhs_h.b', 'gshhs_f.b'};
coastline_file = files{Quality};
fid = fopen(coastline_file, 'r', 'b');
% read in all of the data blocks from the file
tot_points = 1;
lons = NaN;
lats = NaN;
while true
    % read the header of the block
    [polygon_id, count] = fread(fid, 1, 'int32'); % Unique polygon id number, starting at 0
    if (count < 1)
        break;
    end
    num_points = fread(fid, 1, 'int32');  % Number of points in this polygon
    flag = fread(fid, 1, 'int32');
    level = bitand(flag, 255); % 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake
    version = bitand(bitshift(flag,-8),255); % Should be 7 for GSHHS release 7 (i.e., version 2.0)
    greenwich = bitget(flag,17); % Greenwich is 1 if Greenwich is crossed
    source = bitget(flag,25); % 0 = CIA WDBII, 1 = WVS
    river = bitget(flag,26); % 0 = not set, 1 = river-lake and level = 2
    % Min/max extent stored in file in micro-degrees
    west  = fread(fid, 1, 'int32') ./ 1e6; % degrees
    east  = fread(fid, 1, 'int32') ./ 1e6; % degrees
    south = fread(fid, 1, 'int32') ./ 1e6; % degrees
    north = fread(fid, 1, 'int32') ./ 1e6; % degrees
    area  = fread(fid, 1, 'int32') ./10;   % Area of polygon in km^2
    area_full  = fread(fid, 1, 'int32') ./10; % Area of original full-resolution polygon in km^2
    container = fread(fid, 1, 'int32'); % Id of container polygon that encloses this polygon (-1 if none)
    % Id of ancestor polygon in the full resolution set that was the source of this polygon (-1 if none)
    ancestor = fread(fid, 1, 'int32');
    % make sure that this is a coast for land of a reasonable size
    if (level <= maxLevel) && (area(1) > minArea)
        if checkAOI(east,west,north,south)
            % read in the lons/lats of all the data points
            data = fread(fid, [2,num_points] , 'int32') ./ 1e6;

            % shift the range of the longitudes
            hi = find(data(1,:) > 180);
            data(1,hi) = data(1,hi) - 360;

            % check for wrap around
            wraps = [find(abs(diff(data(1, :))) > 180), num_points];
            num_wraps = length(wraps);
            unwrap = data(:,1:wraps(1));
            if num_wraps > 1
                check = zeros(num_wraps-1,1);
                for index = 2:num_wraps
                    check(index-1) = checkAOI(max(data(2,wraps(index-1)+1:wraps(index))),...
                        min(data(2,wraps(index-1)+1:wraps(index))),...
                        max(data(1,wraps(index-1)+1:wraps(index))),...
                        min(data(1,wraps(index-1)+1:wraps(index))));

                    unwrap = [unwrap, [NaN;NaN], data(:,wraps(index-1)+1:wraps(index))];
                end
                if ~any(check)
                    continue
                end
            end
            data = unwrap;
            num_points = num_points + num_wraps - 1;

            % concatenate them to the rest of the list
            tot_points = tot_points + num_points + 1;
            lons = [lons, data(1, :), NaN];
            lats = [lats, data(2, :), NaN];
        else
            % skip over the data points
            fseek(fid, 8.*num_points, 'cof');
        end
    else
        % skip over the data points
        fseek(fid, 8.*num_points, 'cof');
    end
end

% close the data file
fclose(fid);

    function check = checkAOI(east,west,north,south)
        if LatLonCheck
            if west > 180
                west = west-360;
            end
            if east > 180
                east = east-360;
            end
            %longitudinal exclusion
            LonExclude = (west < LatLon(2) && east < LatLon(2)) || (west > LatLon(4) && east > LatLon(4));
            %latitudinal exclusion
            LatExclude = (south < LatLon(1) && north < LatLon(1)) || (south > LatLon(3) && north > LatLon(3));
            if LonExclude || LatExclude
                check = 0;
            else
                %longitudinal crossing
                LonCross = (west < LatLon(2) && east > LatLon(2)) || (west < LatLon(4) && east > LatLon(4));
                %latitudinal crossing
                LatCross = (south < LatLon(1) && north > LatLon(1)) || (south < LatLon(3) && north > LatLon(3));
                if LonCross && ~LonExclude
                    check = 1;
                else if LatCross && ~LatExclude
                        check = 1;
                    else
                        check = 0;
                    end
                end
            end
        else
            check = 1;
        end
    end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function ni = niftiCheckQto(ni)
% Changes the qto_* or/and sto_* fields in the NIFTI hearder if not properly
% set. Also changes other fields see below.
%
%   ni = niftiCheckQto(ni)
%
% We expect these fields to be set at some coordinates inside the image
% volume. If the NIFTI file was created without setting these files this
% function starts up the fields by setting the center of the coordinate
% system of the image to the middle of the image.
%
% Does a simple sanity-check on the xform. Right now, we check the origin to
% make sure it is well- within the image volume. If not, the qto xform
% fields of the NIFTI struct ni will be adjusted to set the origin to the
% center of the image. The updated struct is returned and a message is
% printed to the command line indicating the fix.
%
% it will update the follwing fields:
%  'qform_code'
%  'qto_xyz'
%  'qto_ijk'
%  'quatern_b'
%  'quatern_c'
%  'quatern_d'
%  'qoffset_x'
%  'qoffset_y'
%  'qoffset_z'
%  'qfac'
%  'pixdim'
%
% (c) Stanford Vista Team 2012

% if(isfield(ni,'data')&&~isempty(niftiGet(ni, 'data')))
%     % sanity-check ni.dim
%     sz = size(niftiGet(ni, 'data'));
%     if(any(ni.dim(1:3)~=sz(1:3)))
%         fprintf('[%s] NIFTI volume dim wrong- setting it to the actual data size.\n',mfilename);
%         ni.dim(1:3) = sz(1:3);
%     end
% end

if(niftiGet(ni, 'qform_code')==0 && niftiGet(ni, 'sform_code')~=0)
    fprintf('[%s] ni.qform_code is zero and sform_code ~=0. Setting ni.qto_* from ni.sto_*...\n',mfilename);
    ni = niftiSetQto(ni, niftiGet(ni, 'sto_xyz'));
end

qto_ijk = niftiGet(ni, 'qto_ijk');
origin = [qto_ijk(1:3,:)*[0 0 0 1]']';

nidim = niftiGet(ni, 'dim');
if(any(origin<2)||any(origin>nidim(1:3)-2))
  [t,r,s,k] = affineDecompose(qto_ijk);
  t = nidim/2;
  fprintf('[%s] NIFTI header origin is at or outside the image volume.\n',mfilename)
  fprintf('[%s] Origin to the image center [%2.3f,%2.3f,%2.3f] pix.\n',mfilename,t(1),t(2),t(3));
  ni = niftiSetQto(ni, inv(affineBuild(t,r,s,k)));
end

return;


text:
function res = read(varargin)
opts.scoresroot = fullfile(hb_path, 'scores_new');
opts = vl_argparse(opts, varargin);
scoresroot = opts.scoresroot;
dirs = utls.listdirs(scoresroot);

res.verification = cell(1, numel(dirs));
res.matching = cell(1, numel(dirs));
res.retrieval = cell(1, numel(dirs));

for resi = 1:numel(dirs)
  res.verification{resi} = readtable(fullfile(scoresroot, dirs{resi}, 'verification.csv'));
  res.matching{resi} = readtable(fullfile(scoresroot, dirs{resi}, 'matching.csv'));
  res.retrieval{resi} = readtable(fullfile(scoresroot, dirs{resi}, 'retrieval.csv'));
end
res.verification = vertcat(res.verification{:});
res.matching = vertcat(res.matching{:});
res.retrieval = vertcat(res.retrieval{:});
end

author: hpatches

repo: hpatches-benchmark


text:
% REREF - convert common reference EEG data to some other common reference
%           or to average reference
% Usage:
%   >> Dataout = reref(data);  % convert all channels to average reference
%   >> [Dataout Chanlocs] = reref(data, refchan, 'key', 'val');
%                              % convert data to new reference with options
% Inputs:
%   data - 2-D or 3-D data matrix (chans,frames*epochs)
%   refchan  - reference channel number(s). There are three possibilities:
%          1) [] -  compute average reference
%          2) [X]: re-reference to channel X
%          2) [X Y Z ...]: re-reference to the average of channel X Y Z ...
%
% Optional inputs:
%   'exclude'    - [integer array] channel indices to exclude from re-referencing
%                  (e.g., event marker channels, etc.)
%   'keepref'    - ['on'|'off'] keep reference channel in output (only usable
%                  when there are several references).
%   'elocs'      - Current data electrode location structure (e.g., EEG.chanlocs).
%   'refloc'     - Reference channel location single element structure or cell array
%                  {'label' theta radius} containing the name and polar coordinates
%                  of the current channel. Including this entry means that
%                  this channel will be included (reconstructed) in the
%                  output.
%
% Outputs:
%   Dataout     - Input data converted to the new reference
%   Chanlocs    - Updated channel locations structure
%
% Notes: 1) The average reference calculation implements two methods
%           (see www.egi.com/Technotes/AverageReference.pdf)
%            V'i = (Vi-Vref) - sum(Vi-Vref)/number_of_electrodes
%
% Authors: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2009-
%          previous version: Arnaud Delorme & Scott Makeig, 1999-2002

% Deprecated inputs:
% These inputs are still accepted but not processed. The function returns
% accurate results irrespective of the entry of these options.
%   'refstate  ' - ['common'|'averef'|[indices]] Current reference condition,
%                  ('averef') = average reference; ('common' or 0) = common
%                  reference. [indices] designate the current reference channel
%                  or channels if present in the data {default: 'common'}
%   'method'     - ['standard'|'withref'] Do not ('standard') or do ('withref')
%                  include reference channel data in output {def: 'standard'}.
%                  Note: Option 'withref' not possible when multiple ref channel
%                  indices are given as argument to 'refstate' (below).
%
% ICA inputs:
% These inputs are still accepted but not the ICA conversion is now
% performed from within POP_REREF
%   'icaweights' - ICA weight matrix. Note: If this is ICA weights*sphere,
%                  then the 'icasphere' input below should be [] or identity.
%   'icasphere'  - ICA sphere matrix (if any)
%   'icachansind' - Indices of the channels used in ICA decomposition
%
% Outputs:
%   Wout        - ICA weight matrix (former icaweights*icasphere)
%                 converted to new data reference
%   Sout        - ICA sphere matrix converted to an identity matrix
%   ICAinds     - New indices of channels used in ICA decomposition
%   meandata    - (1,dataframes) means removed from each data point
%
%        2) In conversion of the weight matrix to a new reference
%           where WS = Wts*Sph  and  ica_act = WS*data, then
%              data = inv(WS)*ica_act;
%           If R*data are the re-referenced data,
%              R*data= R*inv(WS)*ica_act;
%           And Wout = inv(R*inv(WS));
%           Now, Sout = eye(length(ICAinds));
%           The re-referenced ICA component maps are now the
%           columns of inv(Wout), and the icasphere matrix, Sout,
%           is an identity matrix. Note: INV -> PINV when
%           PCA dimension reduction is used during ICA decomposition.

% Copyright (C) 1999 Scott Makeig, SCCN/INC/UCSD, scott@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.

% 12/16/99 Corrected denomiator on the suggestion of Ian Nimmo-Smith, Cambridge UK
% 01-25-02 reformated help & license -ad

function [data, Elocs, morechans, W, S, icachansind, meandata] = reref(data, ref, varargin)

if nargin<1
  help reref
  return
end
if nargin < 2
    ref = [];
end

% check inputs
% ------------
g = finputcheck(varargin, { 'icaweight'   'real'    []          [];
                            'icaweights'  'real'    []          [];
                            'icasphere'   'real'    []          [];
                            'icachansind' 'integer'    []       [];
                            'interpchan'  {''}      []          [];
                            'method'     'string'  { 'standard','withref' }  'standard';
                            'refstate'   { 'string','integer' } { { 'common','averef' } [1 size(data,1)] }     'common'; % ot used but kept for backward compatib.
                            'exclude'    'integer' [1 size(data,1)]          [];
                            'refloc'     { 'cell','struct' }  { [] [] }   {};
                            'keepref'    'string'  {'on','off' }             'off';
                            'elocs'      {'integer','struct'}  []            [] });
if ischar(g), error(g); end
if ~isempty(g.icaweight)
    g.icaweights = g.icaweight;
end
if ~isempty(g.icaweights)
    if isempty(g.icachansind),
        g.icachansind = [1:size(g.icaweights,2)];
        disp('Warning: reref() output has changed slightly since EEGLAB 5.02');
        disp('         the 4th output argument is the indices of channels used for ICA instead');
        disp('         of the mean reference value (which is now output argument 5)');
    end
end
if ischar(ref), ref = { ref }; end
if iscell(ref), ref = eeg_chaninds(g.elocs, ref); end
if ~isempty(ref)
    if ref > size(data,1)
        error('reference channel index out of range');
    end
end

[dim1, dim2, dim3] = size(data);
data = reshape(data, dim1, dim2*dim3);

% single reference not present in the data
% add it as blank data channel at the end
% ----------------------------------------
if ~isempty(g.refloc) == 1
    if ~isempty(g.elocs)
        if iscell(g.refloc)
            data(end+1,:) = 0;
            g.elocs(end+1).labels = g.refloc{1};
            g.elocs(end  ).theta  = g.refloc{2};
            g.elocs(end  ).radius = g.refloc{3};
        else
            data(end+length(g.refloc),:) = 0;
            for iLocs = 1:length(g.refloc)
                g.elocs(end+1).labels = g.refloc(iLocs).labels;
                fieldloc = fieldnames(g.elocs);
                for ind = 1:length(fieldloc)
                    g.elocs(end) = setfield(g.elocs(end), fieldloc{ind}, getfield(g.refloc(iLocs), fieldloc{ind}));
                end
            end
        end
    end
    [dim1 dim2 dim3] = size(data);
end

% exclude some channels
% ---------------------
chansin   = setdiff_bc([1:dim1], g.exclude);
nchansin  = length(chansin);

% return mean data
% ----------------
if nargout > 4
    meandata = sum(data(chansin,2))/nchansin;
end

% generate rereferencing matrix
% -----------------------------
if 0 % alternate code - should work exactly the same
    if isempty(ref)
        ref=chansin; % average reference
    end % if
    chansout=chansin;
    data(chansout,:)=data(chansout,:)-ones(nchansin,1)*mean(data(ref,:),1);
else
    if ~isempty(ref) % not average reference
        refmatrix = eye(nchansin); % begin with identity matrix
        tmpref = ref;
        for index = length(g.exclude):-1:1
            tmpref(find(g.exclude(index) < tmpref)) = tmpref(find(g.exclude(index) < tmpref))-1;
        end
        for index = 1:length(tmpref)
            refmatrix(:,tmpref(index)) = refmatrix(:,tmpref(index))-1/length(tmpref);
        end
    else % compute average reference
        refmatrix = eye(nchansin)-ones(nchansin)*1/nchansin;
    end
    chansout = chansin;
    data(chansout,:) = refmatrix*data(chansin,:);
end

% change reference in elocs structure
% -----------------------------------
if ~isempty(g.elocs)
    if isempty(ref)
        for ind = chansin
            g.elocs(ind).ref = 'average';
        end
    else
        reftxt = { g.elocs(ref).labels };
        if length(reftxt) == 1, reftxt = reftxt{1};
        else
            reftxt = cellfun(@(x)([x ' ']), reftxt, 'uniformoutput', false);
            reftxt = [ reftxt{:} ];
        end
        for ind = chansin
            g.elocs(ind).ref = reftxt;
        end
    end
end

% remove reference
% ----------------
morechans = [];
if strcmpi(g.keepref, 'off')
    data(ref,:) = [];
    if ~isempty(g.elocs)
        morechans = g.elocs(ref);
        g.elocs(ref) = [];
    end
end

data = reshape(data, size(data,1), dim2, dim3);

% treat optional ica parameters
% -----------------------------
W = []; S = []; icachansind = [];
if ~isempty(g.icaweights)
    disp('Warning: This function does not process ICA array anymore, use the pop_reref function instead');
end
Elocs = g.elocs;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [] = showBadChannelDifferences(noisy1, noisy2)
%% Output the differences in bad channels between two runs
    numberDatasets = length(noisy1);
    for k = 1:numberDatasets
        badChans1 = noisy1(k).badChannelNumbers;
        badChans2 = noisy2(k).badChannelNumbers;
        u = union(badChans1, badChans2);
        inter = intersect(badChans1, badChans2);
        if length(u) ~= length(inter)
            theDiff = setdiff(u, inter);
            fprintf('%d: [%s]\n', k, num2str(theDiff(:)'));
        end
    end

author: VisLab

repo: EEG-Clean-Tools


text:
function [J, grad] = costFunctionReg(theta, X, y, lambda)
%COSTFUNCTIONREG Compute cost and gradient for logistic regression with regularization
%   J = COSTFUNCTIONREG(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta

H_theta = sigmoid(X * theta);

J = -1/m * sum(y .* log(H_theta) + (1 - y) .* log(1 - H_theta)) + ...
    lambda/(2 * m) * norm(theta([2:end]))^2;

% extra term for gradient
G = lambda/m .* theta;
G(1) = 0;

grad = 1/m .* X' * (H_theta - y) + G;

% =============================================================

end

author: rmarquis

repo: coursera-machinelearning


text:
function [ Z, H, dnorm ] = deep_seminmf ( X, layers, varargin )

% Process optional arguments
pnames = { ...
    'z0' 'h0' 'bUpdateH' 'bUpdateLastH' 'maxiter' 'TolFun', 'verbose', 'bUpdateZ', 'cache' ...
};

num_of_layers = numel(layers);

Z = cell(1, num_of_layers);
H = cell(1, num_of_layers);

dflts  = {0, 0, 1, 1, 500, 1e-5, 1, 1, 1};

[z0, h0, bUpdateH, bUpdateLastH, maxiter, tolfun, verbose, bUpdateZ, cache] = ...
        internal.stats.parseArgs(pnames,dflts,varargin{:});

if  ~iscell(h0)
    for i_layer = 1:length(layers)
        if i_layer == 1
            % For the first layer we go linear from X to Z*H, so we use id
            V = X;
        else
            V = H{i_layer-1};
        end

        if verbose
            display(sprintf('Initialising Layer #%d with k=%d with size(V)=%s...', i_layer, layers(i_layer), mat2str(size(V))));
        end

        % For the later layers we use nonlinearities as we go from
        % g(H_{k-1}) to Z*H_k
        [Z{i_layer}, H{i_layer}, init_err{i_layer}] = ...
             seminmf(V, ...
                 layers(i_layer), ...
                 'maxiter', maxiter, ...
                 'bUpdateH', true, 'bUpdateZ', bUpdateZ, 'verbose', verbose);
                 %'bUpdateH', true, 'bUpdateZ', bUpdateZ, 'verbose', verbose, 'save', cache);
    end

else
    Z=z0;
    H=h0;

    if verbose
        display('Skipping initialization, using provided init matrices...');
    end
end

dnorm0 = cost_function(X, Z, H);
dnorm = dnorm0;

if verbose
    display(sprintf('#%d error: %f', 0, dnorm0));
end

%%% Error Propagation
if verbose
    display('Finetuning...');
end
H_err = cell(1, num_of_layers);

for iter = 1:maxiter
    H_err{numel(layers)} = H{numel(layers)};
    for i_layer = numel(layers)-1:-1:1
        H_err{i_layer} = Z{i_layer+1} * H_err{i_layer+1};
    end

    for i = 1:numel(layers)
        if bUpdateZ
            try
                if i == 1
                    Z{i} = X  * pinv(H_err{1});
                else
                    Z{i} = pinv(D') * X * pinv(H_err{i});
                end
            catch
                display(sprintf('Convergance error %f. min Z{i}: %f. max %f', norm(Z{i}, 'fro'), min(min(Z{i})), max(max(Z{i}))));
            end
        end

        if i == 1
            D = Z{1}';
        else
            D = Z{i}' * D;
        end

        if bUpdateH && (i < numel(layers) || (i == numel(layers) && bUpdateLastH))
            A = D * X;
            Ap = (abs(A)+A)./2;
            An = (abs(A)-A)./2;

            B = D * D';

            Bp = (abs(B)+B)./2;
            Bn = (abs(B)-B)./2;

            H{i} = H{i} .* sqrt((Ap + Bn * H{i}) ./ max(An + Bp * H{i}, 1e-10));
        end
    end

    assert(i == numel(layers));

    dnorm = cost_function(X, Z, H);

    if verbose
        display(sprintf('#%d error: %f', iter, dnorm));
    end

    assert(dnorm <= dnorm0 + 1, ...
        sprintf('Rec. error increasing! From %f to %f. (%d)', ...
        dnorm0, dnorm, iter) ...
    );

    if dnorm0-dnorm <= tolfun*max(1,dnorm0)
        if verbose
            display( ...
                sprintf('Stopped at %d: dnorm: %f, dnorm0: %f', ...
                    iter, dnorm, dnorm0 ...
                ) ...
            );
        end
        break;
    end

    dnorm0 = dnorm;
end
end

function error = cost_function(X, Z, H)
    error = norm(X - reconstruction(Z, H), 'fro');
end

function [ out ] = reconstruction( Z, H )

    out = H{numel(H)};

    for k = numel(H) : -1 : 1;
        out =  Z{k} * out;
    end

end

author: andrewssobral

repo: lrslibrary


text:
% The COBRAToolbox: testDeletionStudy.m
%
% Purpose:
%     - tests the basic functionality of singleGeneDeletion/doubleGeneDeletion/singleRxnDeletion
%       Makes sure that these resultant values are correct:
%       singleGeneDeletion - hasEffect, delRxns, grRateKO, grRateWT
%       doubleGeneDeletion - grRatioDble, grRateKO, grRateWT
%       singleRXDeletion - hasEffect, delRxn, grRateKO, grRateWT
%       returns 1 if all are correct, else 0
%
% Author:
%     - Original file: Joseph Kang 11/16/09
%     - CI integration: Laurent Heirendt

global CBTDIR

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testDeletionStudy'));
cd(fileDir);

tol = 1e-6;

%load model
model = getDistributedModel('ecoli_core_model.mat');

% list of solver packages
solverPkgs = {'tomlab_cplex', 'gurobi', 'glpk'};
%Load reference data.
load('rxnDeletionData.mat');

for k = 1:length(solverPkgs)

    fprintf(' -- Running testfindBlockedReaction using the solver interface: %s ... ', solverPkgs{k});

    solverLPOK = changeCobraSolver(solverPkgs{k}, 'LP', 0);

    if solverLPOK

        fprintf('\n*** Test basic single gene deletion: ***\n\n');

        %deleting gene for 'ENO')
        [grRatio, grRateKO, grRateWT, hasEffect, delRxns, fluxSolution] = singleGeneDeletion(model, 'FBA', {model.genes{1:4},'b2779'});

        % check if correct hasEffect value
        assert(isequal(hasEffect,hasEffectSD))

        % check if correctly deleted reactions
        assert(isequal(delRxns,delRxnsSD))

        % check if correct grRateKO value
        assert(all(abs(grRateKO-grRateSdKO) < tol))

        % check if correct grRateWT value
        assert(abs(grRateWT-grRateWTRef) < tol)

        %Now, we combine gene 1 and two. 1 has no effect, so 1 and 2 should
        %yield the same as 2.
        modelForUTest = model;
        modelForUTest.genes([1,5]) = strcat(model.genes(1),{'.1','.2'});
        targetValues = [1 2 3 4 ];
        %Check functionality of uniqueGene Flag
        [grRatio, grRateKO, grRateWT, hasEffect, delRxns] = singleGeneDeletion(modelForUTest,'FBA',modelForUTest.genes([1 2 3 4]),true,true);
                % check if correct hasEffect value
        assert(isequal(hasEffect,hasEffectSD(targetValues)))

        % check if correctly deleted reactions
        assert(isequal(delRxns,delRxnsSD(targetValues)))

        % check if correct grRateKO value
        assert(all(abs(grRateKO-grRateSdKO(targetValues)) < tol))

        % check if correct grRateWT value
        assert(abs(grRateWT-grRateWTRef) < tol)

        [grRatioDble, grRateKO, grRateWT] = doubleGeneDeletion(model, 'FBA', model.genes(1:4), {'b2779','b2287'});
        %Check against reference
        assert(all(all(abs(grRatioDble-grRatioDbRef) < tol)));

        % check if correct grRateDble value
        assert(all(all(abs(grRateKO-grRateDbKO) < tol)));

        % check if correct grRateWT value
        assert(abs(grRateWT - grRateWTRef) < tol);

        %% singleRxnDeletion Test
        fprintf('\n\nStarting singleRxnDeletion test:\n');

        [test_grRatio, test_grRateKO, test_grRateWT, test_hasEffect, test_delRxn]= singleRxnDeletion(model, 'FBA');

        grRatio_Rxn_Ref(isnan(grRatio_Rxn_Ref)) = -1;
        test_grRatio(isnan(test_grRatio)) = -1;

        % check if correct grRatio values
        for i = 1:length(grRatio_Rxn_Ref)
            assert(abs(grRatio_Rxn_Ref(i) - test_grRatio(i)) < tol)
        end

        grRateKO_Rxn_Ref(isnan(grRateKO_Rxn_Ref)) = -1;
        test_grRateKO(isnan(test_grRateKO)) = -1;

        % check if correct grRateKO values
        for i = 1:length(grRateKO_Rxn_Ref)
            assert(abs(grRateKO_Rxn_Ref(i) - test_grRateKO(i)) < tol)
        end

        % check if correct grRateWT values
        assert(abs(grRateWTRef - test_grRateWT) < tol)

        % check if correct delRxn values
        assert(isequal(delRxn_Rxn_Ref, test_delRxn))
    end

    % output a success message
    fprintf('Done.\n');
end

% change back to root folder
cd(currentDir)

author: opencobra

repo: cobratoolbox


% =========================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% Setup NIREP problem (data needs to be obtained
% separately)
%
% =========================================================================

dataset='na06';
example = ['3D-nirep-',dataset];
checkSetupDataFile; if OK, return; end;

% set view options and interpolation options
[viewer,viewOptn] = viewImage('reset','viewImage','imgmontage','colormap','gray(256)','direction','-zyx');

% setup interpolation scheme
imgOptn  = {'imgModel','splineInterMex','regularizer','moments','theta',.01};

% setup transformation used in the parametric part
traOptn  = {'trafo','affine3D'};

% setup distance measure
disOptn  = {'distance','SSD'};

% initialize the regularizer for the non-parametric part
regOptn = {'regularizer','mbHyperElastic','alpha',1,'alphaLength',1,'alphaArea',.1,'alphaVolume',2};

FAIRmessage(mfilename)
load('na01-128x150x128');
load('na01-128x150x128-labels');
load([dataset,'-128x150x128']);
load([dataset,'-128x150x128-labels']);
dataT = double(dataT);
dataR = double(dataR);
dataTl = double(dataTl);
dataRl = double(dataRl);
m = size(dataT);

dataT = 256.*dataT;
dataR = 256.*dataR;
max(dataT(:))

%omega = [0,m(1),0,m(2),0,m(3)];
omega = [0,20,0,23.4375,0,20];
ML = getMultilevel({dataT,dataR},omega,m,'fig',2);
save(outfile,'dataT','dataR','dataTl','dataRl','omega','m','ML');
save(outfile,'-append','viewOptn','imgOptn','traOptn','disOptn','regOptn');
checkSetupDataFile;

% xc       = getCellCenteredGrid(omega,m);
% viewData = @(I) viewImage(imgModel(I,omega,xc),omega,m);

% FAIRfigure(1,'figname',mfilename); clf;
% subplot(1,2,1); viewData(dataT); title('template');
% subplot(1,2,2); viewData(dataR); title('reference');


text:
% ----------------------------------------------------------------------------
% function hfssSaveProject(fid, projectFile, [Overwrite = false])
%
% Description :
% -------------
% This function creates the necessary VB script to save the active HFSS
% project onto a disk file.
%
% Parameters :
% ------------
% fid         - file identifier of the HFSS script file.
% projectFile - full file name of the project file into which the active HFSS
%               project will be saved.
% [Overwrite] - (boolean) if set to true will overwrite any existing file with
%               the same name as that specified by projectFile.
%
%
% Note :
% ------
% The active HFSS project MUST be saved prior to solving it. This is slightly
% disadvantageous because we cannot create and solve projects "on-the-fly" -
% i.e., in memory without creating a project file on the disk. One workaround
% is to write the project to a temporary file and then delete the temporary
% file after the required solution is obtained.
%
% Example :
% ---------
% fid = fopen('myantenna.vbs', 'wt');
% ...
%
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% This file is part of HFSS-MATLAB-API.
%
% HFSS-MATLAB-API is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the Free
% Software Foundation; either version 2 of the License, or (at your option)
% any later version.
%
% HFSS-MATLAB-API 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
% Foobar; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place, Suite 330, Boston, MA  02111-1307  USA
%
% Copyright 2004, Vijay Ramasami (rvc@ku.edu)
% ----------------------------------------------------------------------------

function hfssSaveProject(fid, projectFile, Overwrite)

% arguments processor.
if (nargin < 2)
	error('Insufficient # of arguments !');
elseif (nargin < 3)
	Overwrite = [];
end

% default arguments.
if isempty(Overwrite)
	Overwrite = false;
end

% create the script.
fprintf(fid, '\n');
fprintf(fid, 'oProject.SaveAs _\n');
fprintf(fid, '    "%s", _\n', projectFile);
if (Overwrite)
    fprintf(fid, '    true\n');
else
    fprintf(fid, '    false\n');
end

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function imshow_numeric ( A )

%*****************************************************************************80
%
%% IMSHOW_NUMERIC displays a numeric 2D array as a grayscale image.
%
%  Discussion:
%
%    A numeric 2D array, typically of type "double", cannot be displayed by
%    the Image Processing Toolbox function imshow.
%
%    This is because the data is of the wrong type (signed real numbers
%    rather than unsigned integers) and generally also in an unsuitable
%    range.  Typical image data might range from 0 to 255, for instance.
%
%    This function makes a quick and simple conversion of numeric data to
%    a format of the appropriate type and range for display by imshow,
%    and displays the data as a grayscale image.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(M,N), the numeric array to be displayed.
%
  a_min = min ( min ( A ) );
  a_max = max ( max ( A ) );

  if ( a_min == a_max )
    A = 127;
  else
    A = floor ( 255 * ( A - a_min ) / ( a_max - a_min ) );
  end

  A = uint8 ( A );

  imshow ( A );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function  idsNms = selectBoundingBoxesNonMaxSup( boundingBoxes, scores, varargin )
%selectBoundingBoxesNonMaxSup performs the NMS on the candidate bounding boxes
% boxes should be in [X1,Y1,W,H] format
%
% idsNms = selectBoundingBoxesNonMaxSup( boundingBoxes, scores )
%
% Input:
%   boundingBoxes - double[ numBoxes x 4], each line correponds to the bounding box in [X1,Y1,W,H] format
%   scores - double[numBoxes x 1], scores of the bounding boxes, will be sorted in the decreasing order
%
% Extra parameters:
%   nmsIntersectionOverAreaThreshold - IoA threshold used to select boxes
%   numBoundingBoxMax - maximum number of boxes selected by NMS
%
% Output:
%   idsNms - indices of the bounding boxes selected by NMS

if ~exist('varargin', 'var')
    varargin = {};
end
%% parameters
opts = struct;
opts.numBoundingBoxMax = inf;
opts.nmsIntersectionOverAreaThreshold = 0.3;
opts = vl_argparse(opts, varargin);

% if opts.nmsIntersectionOverAreaThreshold == inf the code will get into the infinite loop
opts.nmsIntersectionOverAreaThreshold = min( opts.nmsIntersectionOverAreaThreshold, 100 );

%% do the job
numBbs = length(scores);
[~, ids] = sort(scores, 'descend');

idsNms = nan( min(opts.numBoundingBoxMax, numel(ids)), 1 );
idsNms(1) = ids(1);

numBbNms = 1;
iBb = 1;
while numBbNms < opts.numBoundingBoxMax && iBb < numBbs
    curIou = inf;
    while max( curIou(:) ) > opts.nmsIntersectionOverAreaThreshold && iBb < numBbs
        iBb = iBb + 1;
        curIou = bbIntersectionOverArea( boundingBoxes( idsNms(1 : numBbNms), : ), boundingBoxes( ids(iBb), : ) );
    end

    if max( curIou(:) ) <= opts.nmsIntersectionOverAreaThreshold
        numBbNms = numBbNms + 1;
        idsNms( numBbNms ) = ids( iBb );
    end
end

idsNms = idsNms(1 : numBbNms);

end

author: aosokin


text:
% minFunc
fprintf('Compiling minFunc files...\n');
mex -outdir minFunc/compiled minFunc/mex/mcholC.c
mex -outdir minFunc/compiled minFunc/mex/lbfgsC.c
mex -outdir minFunc/compiled minFunc/mex/lbfgsAddC.c
% mex -outdir minFunc/compiled CXX='g++' CXXFLAGS='\$CXXFLAGS -O3 -funroll-loops -fopenmp' LDFLAGS="\$LDFLAGS -fopenmp" minFunc/mex/lbfgsProdC.cpp
if ismac()
    mex -v -outdir minFunc/compiled CXX='g++' CXXFLAGS='\$CXXFLAGS -O3' LDFLAGS="\$LDFLAGS" CXXOPTIMFLAGS='-O3' minFunc/mex/lbfgsProdC.cpp
elseif isunix()
    mex -v -outdir minFunc/compiled CXX='g++' CXXFLAGS='\$CXXFLAGS -O3 -fopenmp -funroll-loops' LDFLAGS="\$LDFLAGS -fopenmp" CXXOPTIMFLAGS='-O3' minFunc/mex/lbfgsProdC.cpp
%     mex -v -outdir minFunc/compiled CXX='g++' CXXFLAGS='\$CXXFLAGS -O0 ' LDFLAGS="\$LDFLAGS " CXXOPTIMFLAGS='-O0' minFunc/mex/lbfgsProdC.cpp
%     mex -outdir minFunc/compiled minFunc/mex/lbfgsProdC.cpp
elseif ispc()
    mex -v -outdir minFunc/compiled COMPFLAGS="-DUSE_WINDOWS /openmp /O3 $COMPFLAGS"  minFunc/mex/lbfgsProdC.cpp
end
% mex -v -outdir minFunc/compiled -fopenmp -O3 minFunc/mex/lbfgsProdC.cpp
% mex -v -outdir minFunc/compiled CXX='g++' CXXFLAGS='\$CXXFLAGS -O3 -msse4 -ftree-vectorize -ftree-vectorizer-verbose=5 -ftree-loop-distribution -funroll-all-loops -ftracer -funroll-loops -fopenmp' LDFLAGS="\$LDFLAGS -fopenmp" minFunc/mex/lbfgsProdC.cpp
% mex -outdir minFunc/compiled CXX='g++' CFLAGS='\$CFLAGS -O3 -msse4 -ftree-vectorize -ftree-vectorizer-verbose=5 -ftree-loop-distribution -funroll-all-loops -ftracer -funroll-loops -fopenmp' LDFLAGS="\$LDFLAGS -fopenmp" minFunc/mex/lbfgsProdC.c

%%
% mex -v -outdir minFunc/compiled COMPFLAGS="-DUSE_WINDOWS /openmp /O3 $COMPFLAGS"  minFunc/mex/lbfgsProdC.cpp

repo: pTVreg


text:
function tf=moxunit_util_platform_is_octave
% return true if and only if the platform is GNU/Octave
%
% tf=moxunit_util_platform_is_octave()
%
% Output:
%   tf      true if the platform on which this function is run is
%           GNU/Octave, false otherwise (typically, if the platform is
%           Matlab)
%
% NNO Jan 2014

    persistent cached_tf;

    if islogical(cached_tf)
        tf=cached_tf;
        return;
    end

    tf=logical(exist('OCTAVE_VERSION', 'builtin'));
    cached_tf=tf;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [B, elapse] = BRE_compress(A, model)
%   This is a wrapper function of Binary Reconstructive Embedding testing.
%
%	Usage:
%	[B,elapse] = BRE_compress(A, model)
%
%	      A: Rows of vectors of data points. Each row is sample point
%     model: The model generated by BRE_learn.
%
%	      B: The binary code of the input data A. Each row is sample point
%    elapse: The coding time (testing time).
%
%
%
%   version 2.0 --Nov/2016
%   version 1.0 --Jan/2013
%
%   Written by  Yue Lin (linyue29@gmail.com)
%               Deng Cai (dengcai AT gmail DOT com)
%

tmp_T = tic;

Ktest = A*model.X';
H_query = zeros(size(A,1),size(model.W,2));
for b = 1:size(model.W,2)
    H_query(:,b) = Ktest(:,model.hash_inds(:,b))*model.W(:,b);
end
B = (H_query>0);

elapse = toc(tmp_T);
end

author: ZJULearning

repo: MatlabFunc


text:
classdef SO3FunHandle < SO3Fun
% a class representing a function on the rotation group by an function
% handle
%
% Syntax
%   SO3F = SO3FunHandle(fun)
%
% Input
%  fun - @function_handle
%
% Output
%  SO3F - @SO3FunHandle
%
% Example
%
%   r = orientation.rand;
%   SO3F = SO3FunHandle(@(rot) angle(rot,r))
%
properties
  fun
  % TODO: antipodal wird weder erkannt noch gesetzt/verwendet
  antipodal = false
  SLeft  = specimenSymmetry
  SRight = specimenSymmetry
  bandwidth = 96
end

methods
  function SO3F = SO3FunHandle(fun,varargin)

    SO3F.fun = fun;

    [SRight,SLeft] = extractSym(varargin);
    SO3F.SRight = SRight;
    SO3F.SLeft = SLeft;

  end

  function f = eval(SO3F,rot,varargin)

%     if isa(rot,'orientation')
%       ensureCompatibleSymmetries(SO3F,rot)
%     end

    s = size(rot);
    rot = rot(:);
    f = SO3F.fun(rot);
    if numel(f)==length(rot)
      f = reshape(f,s);
    end
    if isalmostreal(f)
      f = real(f);
    end
  end

end

methods (Static = true)

  SO3F = example(varargin)

end

end

author: mtex-toolbox

repo: mtex

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


text:
function Script_CECLM_of_multi_hyp()

addpath(genpath('../'));

[images, detections, labels] = Collect_menpo_train_imgs('G:\datasets\menpo/');

%% loading the patch experts
[patches, pdm, clmParams] = Load_CECLM_OF();

%% Fitting the model to the provided image
views = [0,0,0; 0,-30,0; 0,30,0; 0,-55,0; 0,55,0; 0,0,30; 0,0,-30; 0,-90,0; 0,90,0; 0,-70,40; 0,70,-40];
views = views * pi/180;
num_views = size(views,1);

clmParams.numPatchIters = 1;

% for recording purposes
experiment.params = clmParams;

num_points = numel(pdm.M)/3;

shapes_all = cell(numel(images), num_views);
labels_all = cell(numel(images), 1);
lhoods = zeros(numel(images),num_views);
all_lmark_lhoods = zeros(num_points, numel(images),num_views);
all_views_used = zeros(numel(images),num_views);
errors_view = zeros(numel(images),num_views);
% Use the multi-hypothesis model, as bounding box tells nothing about
% orientation
verbose = true;
tic

for i=1:numel(images)

    image = imread(images(i).img);
    image_orig = image;

    if(size(image,3) == 3)
        image = rgb2gray(image);
    end

    bbox = squeeze(detections(i,:));

    % Find the best orientation
    for v = 1:size(views,1)
        [shape,~,~,lhood,lmark_lhoods,view_used] = Fitting_from_bb(image, [], bbox, pdm, patches, clmParams, 'orientation', views(v,:));
        shapes_all{i,v} = shape;
        all_views_used(i,v) = view_used;
        all_lmark_lhoods(:, i,v) = lmark_lhoods;
        lhoods(i,v) = lhood;
        errors_view(i,v) = compute_error_menpo_1(labels(i), {shape});
    end

    labels_all{i} = labels{i};

    if(mod(i, 200)==0)
        fprintf('%d done\n', i );
    end

end
toc

experiment.lhoods = lhoods;
experiment.errors_view = errors_view;
experiment.all_views_used = all_views_used;

%%
output_results = ['results/results_ceclm_of.mat'];
save(output_results, 'experiment');
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function sampleScatterMatrix(rxnNames, model, sample, nPoints, fontSize, dispRFlag, rxnNames2)
% Draws a scatterplot matrix with pairwise scatterplots
% for multiple reactions
%
% USAGE:
%
%    sampleScatterMatrix(rxnNames, model, sample, nPoints, fontSize, dispRFlag, rxnNames2)
%
% INPUTS:
%    rxnNames:     Cell array of reaction names to be plotted
%    model:        Model structure
%    sample:       Samples to be analyzed (`nRxns x nSamples`)
%
% OPTIONAL INPUTS:
%    nPoints:      How many sample points to plot (Default `nSamples`)
%    fontSize:     Font size for labels (Default calculated based on
%                  number of reactions)
%    dispRFlag:    Display correlation coefficients (Default false)
%    rxnNames2:    Optional second set of reaction names
%
% EXAMPLE:
%
%    %Plots the scatterplots only between the three reactions listed -
%    %histograms for each reaction will be on the diagonal
%    sampleScatterMatrix({'PFK', 'PYK', 'PGL'}, model, sample);
%    %Plots the scatterplots between each of the first set of reactions and
%    %each of the second set of reactions. No histograms will be shown.
%    sampleScatterMatrix({'PFK', 'PYK', 'PGL'}, model, sample, 100, 10, true, {'ENO','TPI'});
%
% .. Author: - Markus Herrgard 9/14/06

[isInModel, rxnInd] = ismember(rxnNames, model.rxns);
rxnNames = rxnNames(isInModel);
rxnInd = rxnInd(isInModel);
nRxns = length(rxnNames);

%no optional inputs specified
if nargin < 4
    nPoints = size(sample, 2);
    dispRFlag = false;
    nRxns2 = nRxns;
    rxnNames2 = rxnNames;
    nPanels = nRxns * nRxns2;
    fontSize = 10 + ceil(50 / sqrt(nPanels));
    twoSetsFlag = false;
%some optional inputs specified
else
    if isempty(nPoints)
        nPoints = size(sample, 2);
    end
    if isempty(dispRFlag)
        dispRFlag = false;
    end
    if nargin == 7
        [isInModel2, rxnInd2] = ismember(rxnNames2, model.rxns);

        rxnNames2 = rxnNames2(isInModel2);
        rxnInd2 = rxnInd2(isInModel2);

        nRxns2 = length(rxnNames2);
        twoSetsFlag = true;
        nRxns = nRxns + 1;
        nRxns2 = nRxns2 + 1;
    else
        nRxns2 = nRxns;
        rxnNames2 = rxnNames;
        twoSetsFlag = false;
    end
    if isempty(fontSize)
        nPanels = nRxns * nRxns2;
        fontSize = 10 + ceil(50 / sqrt(nPanels));
    end
end
height = 0.8 / nRxns;
width = 0.8 / nRxns2;

clf
showprogress(0, 'Drawing scatterplots ...');
for i = 1:nRxns

    for j = 1:nRxns2
        showprogress(((i - 1) * nRxns+j) / (nRxns * nRxns2));
        left = 0.1 + (j - 1) * width;
        bottom = 0.9 - i * height;
        if twoSetsFlag
            if i == 1 && j == 1
            else
                %subplot(nRxns,nRxns2,(i-1)*nRxns2+j);
                subplot('position', [left bottom width height]);
                if j >1 && i >1
                    sampleScatterPlot(sample, rxnInd2(j - 1), rxnInd(i - 1), nPoints, fontSize, dispRFlag);
                elseif i == 1
                    sampleHistInternal(sample, rxnInd2(j - 1), fontSize);
                elseif j == 1
                    sampleHistInternal(sample, rxnInd(i - 1), fontSize);
                end
                if i == 1
                    xlabel(rxnNames2{j - 1}, 'FontSize', fontSize);
                    set(gca, 'XAxisLocation', 'top');
                end
                if j == 1
                    ylabel(rxnNames{i - 1}, 'FontSize', fontSize);
                end
            end
        else
            if j == i
                %subplot(nRxns,nRxns2,(i-1)*nRxns2+j);
                subplot('position', [left bottom width height]);
                sampleHistInternal(sample, rxnInd(i), fontSize);

            elseif j > i
                %subplot(nRxns,nRxns2,(i-1)*nRxns2+j);
                subplot('position', [left bottom width height]);
                sampleScatterPlot(sample, rxnInd(j), rxnInd(i), nPoints, fontSize, dispRFlag);
            end
            if i == 1
                xlabel(rxnNames2{j}, 'FontSize', fontSize);
                set(gca, 'XAxisLocation', 'top');
            end
            if j == nRxns2
                set(gca, 'YAxisLocation', 'right');
                ylabel(rxnNames{i}, 'FontSize', fontSize);
            end
        end
    end

end

function sampleScatterPlot(sample, id1, id2, nPoints, fontSize, dispRFlag)

selPts = randperm(size(sample, 2));
selPts = selPts(1:nPoints);

plot(sample(id1, selPts),sample(id2, selPts), 'r.');
set(gca, 'YTickLabel', []);
set(gca, 'XTickLabel', []);
maxx = max(sample(id1, :));
maxy = max(sample(id2, :));
minx = min(sample(id1, :));
miny = min(sample(id2, :));
axis([minx maxx miny maxy]);
% Display correlation coefficients
if dispRFlag
    r = corrcoef(sample(id1, :)', sample(id2, :)');
    h = text(minx + 0.66 * (maxx-minx), miny + 0.2 * (maxy - miny), num2str(round(100 * r(1, 2)) / 100));
    set(h, 'FontSize', fontSize - 5);
end

function sampleHistInternal(sample , id, fontSize)

[n, bins] = hist(sample(id , :), 30);
if exist('smooth')
    plot(bins, smooth(bins, n') / sum(n'));
else
    plot(bins, n' / sum(n'));
end
maxx = max(bins);
minx = min(bins);
set(gca, 'XTick', linspace(minx, maxx, 4));
set(gca, 'XTickLabel', round(10 * linspace(minx, maxx, 4)) / 10);
set(gca, 'YTickLabel', []);
set(gca, 'FontSize', fontSize - 5);
axis tight


text:
if and(ssub==1, tsub==1)
    neuron = neuron_full;
    Y = double(data.Y(:, :, sframe+(1:num2read)-1));
    [d1s,d2s, T] = size(Y);
    fprintf('\nThe data has been loaded into RAM. It has %d X %d pixels X %d frames. \nLoading all data requires %.2f GB RAM\n\n', d1s, d2s, T, d1s*d2s*T*8/(2^30));
else
    if exist(nam_mat, 'file')
        [Y, neuron_ds] = neuron_full.load_data(nam_mat, sframe, num2read);
        [d1s,d2s, T] = size(Y);
        fprintf('\nThe data has been downsampled and loaded into RAM. It has %d X %d pixels X %d frames. \nLoading all data requires %.2f GB RAM\n\n', d1s, d2s, T, d1s*d2s*T*8/(2^30));
        neuron = neuron_ds.copy();
    else
        [Y, neuron_ds] = neuron_full.load_data(nam, sframe, num2read);
        [d1s,d2s, T] = size(Y);
        fprintf('\nThe data has been downsampled and loaded into RAM. It has %d X %d pixels X %d frames. \nLoading all data requires %.2f GB RAM\n\n', d1s, d2s, T, d1s*d2s*T*8/(2^30));
        neuron = neuron_ds.copy();
    end
end

author: zhoupc


text:
%compute dtailheadang
function [data,units]=compute_dtailheadang(trx,n)
larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dtailheadang=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    % KB: this angle distance didn't make sense to me, changed to one that
    % made sense to me
    dtailheadang{1,i}=modrange(diff(trx(larva).tailheadang),-pi,pi)./trx(larva).dt;
    % dtailheadang{1,i}=(mod1(trx(larva).tailheadang(2:end)-trx(larva).tailheadang(1:end-1),pi)-pi)./trx(larva).dt;
end

units=parseunits('rad/s');
data=dtailheadang;

author: kristinbranson

repo: JAABA