text:
function sMap = som_read_cod(filename)

%SOM_READ_COD Reads a SOM_PAK format codebook file.
%
% sMap = som_read_cod(filename);
%
%  sMap = som_read_cod('map1.cod');
%
%  Input and output arguments:
%   filename    (string) name of input file
%   sMap        (struct) self-organizing map structure
%
% The file must be in SOM_PAK format. Empty lines and lines starting
% with a '#' are ignored, except the ones starting with '#n'. The strings
% after '#n' are read to field 'comp_names' of the map structure.
%
% For more help, try 'type som_read_cod' or check out online documentation.
% See also SOM_WRITE_COD, SOM_READ_DATA, SOM_WRITE_DATA, SOM_MAP_STRUCT.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_read_cod
%
% PURPOSE
%
% Reads a Self-Organizing Map from an ascii file in SOM_PAK format.
%
% SYNTAX
%
%  sMap = som_read_cod(filename);
%
% DESCRIPTION
%
% This function is offered for compatibility with SOM_PAK, a SOM
% software package in C. It reads map files written in SOM_PAK format.
%
% The SOM_PAK map file format is as follows. The first line must contain
% the input space dimension, lattice type ('rect' or 'hexa'), map grid
% size in x-direction, map grid size in y-direction, and neighborhood
% function ('bubble' or 'gaussian'), in that order. The following lines
% are comment lines, empty lines or data lines.
%
% Each data line contains the weight vector of one map unit and its
% labels. From the beginning of the line, first are values of the vector
% components separated by whitespaces, then labels, again separated by
% whitespaces. The order of map units in the file are one row at a time
% from right to left, from the top to the bottom of the map (x-direction
% first, then y-direction).
%
% Comment lines start with '#'. Comment lines as well as empty lines are
% ignored, except if the comment line starts with '#n'. In that case the
% line should contain names of the vector components separated by
% whitespaces.
%
% In the returned map struct, several fields has to be set to default
% values, since the SOM_PAK file does not contain information on
% them. These include map shape ('sheet'), mask ([1 ... 1]),
% normalizations (none), trainhist (two entries, first with algorithm
% 'init' and the second with 'seq', both with data name 'unknown'),
% possibly also component names ('Var1',...).
%
% REQUIRED INPUT PARAMETERS
%
%  filename   (string) the name of the input file
%
% OUTPUT ARGUMENTS
%
%  sMap       (struct) the resulting SOM struct
%
% EXAMPLES
%
%  sMap = som_read_cod('map1.cod');
%
% SEE ALSO
%
%  som_write_cod    Writes a map struct into a file in SOM_PAK format.
%  som_read_data    Reads data from an ascii file.
%  som_write_data   Writes data struct into a file in SOM_PAK format.
%  som_map_struct   Creates map structs.

% Copyright (c) 1997-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 1.0beta ecco 221097
% Version 2.0beta juuso 151199 250400

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

error(nargchk(1, 1, nargin))  % check no. of input args is correct

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialize variables

lnum           = 0;    % codebook vector counter
comment_start  = '#';  % the char a SOM_PAK command line starts with
comp_name_line = '#n'; % string used to start a special command line,
                       % which contains names of each component

% open input file

fid = fopen(filename);
if fid < 0, error(['Cannot open ' filename]); end

% read header line

ok_cnt = 0;
lin = fgetl(fid); li = lin;
[dim c err n]  = sscanf(li, '%d%[^ \t]'); ok_cnt=ok_cnt+c; li = li(n:end);
[lattice c err n] = sscanf(li,'%s%[^ \t]'); ok_cnt=ok_cnt+c; li = li(n:end);
[msize(2) c err n] = sscanf(li, '%d%[^ \t]'); ok_cnt=ok_cnt+c; li = li(n:end);
[msize(1) c err n] = sscanf(li, '%d%[^ \t]'); ok_cnt=ok_cnt+c; li = li(n:end);
[neigh c err n] = sscanf(li, '%s%[^ \t\n]'); ok_cnt=ok_cnt+c;

if ok_cnt ~= 5
  error([ 'Invalid header line: ' lin ]);
end

% create map struct and set its fields according to header line

munits = prod(msize);
sMap   = som_map_struct(dim, 'msize', msize, ...
			lattice, 'sheet', 'neigh', neigh);
[sT0, ok] = som_set('som_train','algorithm','init','data_name','unknown');
sT1       = som_set('som_train','algorithm','seq','data_name','unknown',...
                    'neigh',neigh,'mask',ones(dim,1));
[sMap, ok, msgs] = som_set(sMap,'name',filename,'trainhist',{sT0,sT1});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% read codebook from the file

codebook   = zeros(munits,dim);
labels     = cell(munits,1);
comp_names = sMap.comp_names;
form       = [repmat('%f',[1 dim-1]) '%f%[^ \t]'];

while 1,
  li = fgetl(fid);                         % read next line
  if ~ischar(li), break, end;               % is this the end of file?

  [data, c, err, n] = sscanf(li, form);
  if c < dim % if there were less numbers than dim on the input file line
    if c == 0
      if strncmp(li, comp_name_line, 2) % component name line?
        li = li(3:end); i = 0; c = 1;
        while c
          [s, c, e, n] = sscanf(li, '%s%[^ \t]');
          if ~isempty(s), i = i + 1; comp_names{i} = s; li = li(n:end); end
        end

        if i ~= dim
          error(['Illegal number of component names: ' num2str(i) ...
                 ' (dimension is ' num2str(dim) ')']);
        end
      elseif ~strncmp(li, comment_start, 1) % not a comment, is it error?
        [s, c, e, n] = sscanf(li, '%s%[^ \t]');
        if c
          error(['Invalid vector on input file line ' ...
                 num2str(lnum+1) ': [' deblank(li) ']']),
        end
      end
    else
      error(['Only ' num2str(c) ' vector components on input file line ' ...
             num2str(lnum+1) ' (dimension is ' num2str(dim) ')']);
    end

  else

    lnum = lnum + 1;                % this was a line containing data vector
    codebook(lnum, 1:dim) = data'; % add data to struct

    % read labels

    if n < length(li)
      li = li(n:end);
      i = 0; n = 1; c = 1;
      while c
        [s, c, e, n_new] = sscanf(li(n:end), '%s%[^ \t]');
        if c, i = i + 1; labels{lnum, i} = s; n = n + n_new - 1; end
      end
    end
  end
end

% close the input file

if fclose(fid) < 0
  error(['Cannot close file ' filename]);
else
  fprintf(2, '\rmap read ok         \n');
end

% check that the number of lines read was correct

if lnum ~= munits
  error(['Illegal number of map units: ' num2str(lnum) ' (should be ' num2str(munits) ').']);
end

% set values

% in SOM_PAK the xy-indexing is used, while in Matlab ij-indexing
% therefore, the codebook vectors have to be reorganized

order = reshape([1:munits],msize);
order = reshape(order',[munits 1]);
codebook(order,:) = codebook;
labels(order,:) = labels;

sMap.codebook   = codebook;
sMap.labels     = labels;
sMap.comp_names = comp_names;

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: ilarinieminen

repo: SOM-Toolbox


text:
%{
    Tests the "analysis" form of basis pursuit de-noising

    min_x ||Wx||_1
s.t.
    || A(x) - b || <= eps

The solvers solve a regularized version, using
    ||Wx||_1 + mu/2*||x-x_0||_2^2

see also test_sBP.m and test_sBPDN.m

%}

% Before running this, please add the TFOCS base directory to your path
myAwgn = @(x,snr) x + ...
        10^( (10*log10(sum(abs(x(:)).^2)/length(x(:))) - snr)/20 )*randn(size(x));

% Try to load the problem from disk
fileName = fullfile('reference_solutions','basispursuit_W_problem1_smoothed_noisy');
randn('state',34324);
rand('state',34324);

N = 512;
M = round(N/2);
K = round(M/5);

A = randn(M,N);
x = zeros(N,1);

% introduce a sparsifying transform "W"
d  = round(4*N);     % redundant
Wf = @(x) dct(x,d);  % zero-padded DCT
downsample  = @(x) x(1:N,:);
Wt = @(y) downsample( idct(y) );    % transpose of W
W  = Wf(eye(N));     % make an explicit matrix for CVX
if exist([fileName,'.mat'],'file')
    load(fileName);
    fprintf('Loaded problem from %s\n', fileName );

else

    % Generate a new problem

    % the signal x consists of several pure tones at random frequencies
    for k = 1:K
        x = x + randn()*sin( rand()*pi*(1:N) + 2*pi*rand() ).';
    end

    b_original = A*x;
    snr = 40;  % SNR in dB
    b   = myAwgn(b_original,snr);
    EPS = norm(b-b_original);
    x_original = x;

    mu = .01*norm(Wf(x),Inf);
    x0 = zeros(N,1);

    % get reference via CVX
    tic
    cvx_begin
        cvx_precision best
        variable xcvx(N,1)
        minimize norm(W*xcvx,1) + mu/2*sum_square(xcvx-x0)
        subject to
            norm(A*xcvx - b ) <= EPS
    cvx_end
    time_IPM = toc;
    x_ref = xcvx;
    obj_ref = norm(W*x_ref,1) + mu/2*sum_square(x_ref-x0);

    save(fileName,'x_ref','b','x_original','mu',...
        'EPS','b_original','obj_ref','x0','time_IPM','snr','d');
    fprintf('Saved data to file %s\n', fileName);

end

[M,N]           = size(A);
norm_x_ref      = norm(x_ref);
norm_x_orig     = norm(x_original);
er_ref          = @(x) norm(x-x_ref)/norm_x_ref;
er_signal       = @(x) norm(x-x_original)/norm_x_orig;
resid           = @(x) norm(A*x-b)/norm(b);  % change if b is noisy

fprintf('\tA is %d x %d\n', M, N );
fprintf('\tl1-norm solution and original signal differ by %.2e (mu = %.2e)\n', ...
    norm(x_ref - x_original)/norm(x_original),mu );

%% Call the TFOCS solver
objF            = @(x)norm(W*x,1) + mu/2*norm(x-x0).^2;
infeasF         = @(x)norm(A*x-b) - EPS;
er              = er_ref;  % error with reference solution (from IPM)
opts            = [];
opts.errFcn     = { @(f,dual,primal) er(primal), ...
                    @(f,dual,primal) obj_ref - f, ...
                    @(f,dual,primal) infeasF(primal) };
opts.maxIts     = 1000;
opts.tol        = 1e-10;
% opts.normA2     = norm(A*A');
% opts.normW2     = norm(W'*W);
z0  = [];   % we don't have a good guess for the dual
tic;
[ x, out, optsOut ] = solver_sBPDN_W( A, W, b, EPS, mu, x0, z0, opts );
time_TFOCS = toc;
fprintf('x is sub-optimal by %.2e, and infeasible by %.2e\n',...
    objF(x) - obj_ref, infeasF(x) );

fprintf('Solution has %d nonzeros.  Error vs. IPM solution is %.2e\n',...
    nnz(x), er(x) );

% Check that we are within allowable bounds
if out.err(end,1) < 1e-4
    disp('Everything is working');
else
    error('Failed the test');
end

%% Doing it "by hand" (for debugging), just for 200 iterations
objF            = @(x)norm(W*x,1) + mu/2*norm(x-x0).^2;
infeasF         = @(x)norm(A*x-b) - EPS;
er              = er_ref;  % error with reference solution (from IPM)
% er              = er_signal;
opts            = [];
opts.errFcn     = { @(f,dual,primal) er(primal), ...
                    @(f,dual,primal) obj_ref - f, ...
                    @(f,dual,primal) infeasF(primal) };
opts.maxIts     = 200;
opts.printEvery = 10;
opts.tol        = 1e-10;
z0      = [];
proxScale = sqrt(norm(W'*W)/norm(A*A'));
scale   = 1;
prox    = { prox_l2( EPS/scale ), proj_linf(proxScale) };
affineF = {A/scale,-b/scale;W/proxScale,0};
[ x, out, optsOut ] = tfocs_SCD( [], affineF, prox, mu, x0, z0, opts );

% Check that we are within allowable bounds
if out.err(end,1) < 1e-2
    disp('Everything is working');
else
    error('Failed the test');
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [mri] = ft_read_mri(filename)

% FT_READ_MRI reads anatomical and functional MRI data from different
% file formats. The output data is structured in such a way that it is
% comparable to a FieldTrip source reconstruction.
%
% Use as
%   [mri] = ft_read_mri(filename)
%
% The output MRI may have a homogenous transformation matrix that converts
% the coordinates of each voxel (in xgrid/ygrid/zgrid) into head
% coordinates.
%
% See also FT_READ_DATA, FT_READ_HEADER, FT_READ_EVENT

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

% test for the presence of some external functions from other toolboxes
hasmri  = ft_hastoolbox('mri');     % from Darren Weber, see http://eeg.sourceforge.net/
hasspm2 = ft_hastoolbox('SPM2');    % see http://www.fil.ion.ucl.ac.uk/spm/
hasspm5 = ft_hastoolbox('SPM5');    % see http://www.fil.ion.ucl.ac.uk/spm/
hasspm8 = ft_hastoolbox('SPM8');    % see http://www.fil.ion.ucl.ac.uk/spm/
hasspm = (hasspm2 || hasspm5 || hasspm8);
hasafni = ft_hastoolbox('afni');    % see http://afni.nimh.nih.gov/

% test whether the file exists
if ~exist(filename)
  error(sprintf('file ''%s'' does not exist', filename));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ft_filetype(filename, 'ctf_mri')
  [img, hdr] = read_ctf_mri(filename);
  transform = hdr.transformMRI2Head;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ft_filetype(filename, 'ctf_mri4')
  [img, hdr] = read_ctf_mri4(filename);
  transform = hdr.transformMRI2Head;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ft_filetype(filename, 'ctf_svl')
  [img, hdr] = read_ctf_svl(filename);
  transform = hdr.transform;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ft_filetype(filename, 'asa_mri')
  [img, seg, hdr] = read_asa_mri(filename);
  transform = hdr.transformMRI2Head;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ft_filetype(filename, 'minc')
  if ~(hasspm2 || hasspm5)
    error('the SPM2 or SPM5 toolbox is required to read *.mnc files');
  end
  % use the functions from SPM
  hdr = spm_vol_minc(filename);
  img = spm_read_vols(hdr);
  transform = hdr.mat;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ft_filetype(filename, 'nifti')
  if ~(hasspm5 || hasspm8)
    error('the SPM5 or SPM8 toolbox is required to read *.nii files');
  end
  % use the functions from SPM
  hdr = spm_vol_nifti(filename);
  img = spm_read_vols(hdr);
  transform = hdr.mat;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif (ft_filetype(filename, 'analyze_img') || ft_filetype(filename, 'analyze_hdr')) && hasspm
  % use the image file instead of the header
  filename((end-2):end) = 'img';
  % use the functions from SPM to read the Analyze MRI
  hdr = spm_vol(filename);
  img = spm_read_vols(hdr);
  transform = hdr.mat;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif (ft_filetype(filename, 'analyze_hdr') || ft_filetype(filename, 'analyze_img')) && hasmri
  % use the functions from Darren Weber's mri_toolbox to read the Analyze MRI
  avw = avw_img_read(filename, 0); % returned volume is LAS*
  img = avw.img;
  hdr = avw.hdr;
  % The default Analyze orientation is axial unflipped (LAS*), which means
  % that the resulting volume is according to the radiological convention.
  % Most other fMRI and EEG/MEG software (except Mayo/Analyze) uses
  % neurological conventions and a right-handed coordinate system, hence
  % the first axis of the 3D volume (right-left) should be flipped to make
  % the coordinate system comparable to SPM
  warning('flipping 1st dimension (L-R) to obtain volume in neurological convention');
  img = flipdim(img, 1);
  % FIXME: here I should also implement a homogenous transformation matrix,
  % using the voxel dimensions that are specified in hdr.dime.pixdim

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif (ft_filetype(filename, 'afni_brik') || ft_filetype(filename, 'afni_head')) && hasafni
  [err, img, hdr, ErrMessage] = BrikLoad(filename);
  if err
    error('could not read AFNI file');
  end

  % FIXME: this should be checked, but I only have a single BRIK file
  % construct the homogenous transformation matrix that defines the axes
  warning('homogenous transformation might be incorrect for AFNI file');
  transform        = eye(4);
  transform(1:3,4) = hdr.ORIGIN(:);
  transform(1,1)   = hdr.DELTA(1);
  transform(2,2)   = hdr.DELTA(2);
  transform(3,3)   = hdr.DELTA(3);

  % FIXME: I am not sure about the "RAI" image orientation
  img = flipdim(img,1);
  img = flipdim(img,2);
  dim = size(img);
  transform(1,4) = -dim(1) - transform(1,4);
  transform(2,4) = -dim(2) - transform(2,4);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ft_filetype(filename, 'neuromag_fif') && ft_hastoolbox('mne')
  % use the mne functions to read the Neuromag MRI
  hdr = fiff_read_mri(filename);
  img = cat(3, hdr.slices.data);
  hdr.slices = rmfield(hdr.slices, 'data'); % remove the image data to save memory
  % hmm, which transformation matrix should I use?
  if isfield(hdr, 'voxel_trans') && issubfield(hdr.voxel_trans, 'trans')
    transform = hdr.voxel_trans.trans;
  elseif isfield(hdr, 'trans') && issubfield(hdr.trans, 'trans')
    transform = hdr.trans.trans;
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ft_filetype(filename, 'neuromag_fif') && ft_hastoolbox('meg_pd')
  % use the meg_pd functions to read the Neuromag MRI
  [img,coords] = loadmri(filename);
  dev = loadtrans(filename,'MRI','HEAD');
  transform  = dev*coords;
  hdr.coords = coords;
  hdr.dev    = dev;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ft_filetype(filename, 'neuromag_fif')
  error('reading MRI data from a fif file requires either the MNE toolbox or the meg_pd toolbox to be installed');

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ft_filetype(filename, 'dicom')
  % this uses the Image processing toolbox
  % the DICOM file probably represents a stack of slices, possibly even multiple volumes
  orig = dicominfo(filename);
  dim(1) = orig.Rows;
  dim(2) = orig.Columns;

  [p, f] = fileparts(filename);

  % this works for the Siemens scanners at the FCDC
  tok = tokenize(f, '.');
  for i=5:length(tok)
    tok{i} = '*';
  end
  filename = sprintf('%s.', tok{:});  % reconstruct the filename with wildcards and '.' between the segments
  filename = filename(1:end-1);       % remove the last '.'
  dirlist  = dir(fullfile(p, filename));
  dirlist  = {dirlist.name};

  if length(dirlist)==1
    % try something else to get a list of all the slices
    dirlist = dir(fullfile(p, '*'));
    dirlist = {dirlist(~[dirlist.isdir]).name};
  end

  keep = false(1, length(dirlist));
  for i=1:length(dirlist)
    filename = char(fullfile(p, dirlist{i}));
    if ~ft_filetype(filename, 'dicom')
      keep(i) = false;
      fprintf('skipping ''%s'' because of incorrect filetype\n', filename);
    end
    % read the header information
    info     = dicominfo(filename);
    if info.SeriesNumber~=orig.SeriesNumber
      keep(i) = false;
      fprintf('skipping ''%s'' because of different SeriesNumber\n', filename);
    else
      keep(i) = true;
      hdr(i)  = info;
    end
  end
  % remove the files that were skipped
  hdr     = hdr(keep);
  dirlist = dirlist(keep);

  % pre-allocate enough space for the subsequent slices
  dim(3) = length(dirlist);
  img    = zeros(dim(1), dim(2), dim(3));
  for i=1:length(dirlist)
    filename = char(fullfile(p, dirlist{i}));
    fprintf('reading image data from ''%s''\n', filename);
    img(:,:,i) = dicomread(hdr(i));
  end

  % reorder the slices
  [z, indx]   = sort(cell2mat({hdr.SliceLocation}));
  hdr = hdr(indx);
  img = img(:,:,indx);

  try
    % construct a homgeneous transformation matrix that performs the scaling from voxels to mm
    dx = hdr(1).PixelSpacing(1);
    dy = hdr(1).PixelSpacing(2);
    dz = hdr(2).SliceLocation - hdr(1).SliceLocation;
    transform = eye(4);
    transform(1,1) = dx;
    transform(2,2) = dy;
    transform(3,3) = dz;
  end
elseif ft_filetype(filename, 'freesurfer_mgz')
  ft_hastoolbox('freesurfer', 1);
  tmp = MRIread(filename);
  img = permute(tmp.vol, [2 1 3]); %FIXME although this is probably correct
  %see the help of MRIread, anecdotally columns and rows seem to need a swap
  %in order to match the transform matrix (alternatively a row switch of the
  %latter can be done)
  hdr = rmfield(tmp, 'vol');
  transform = tmp.vox2ras1;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
  error(sprintf('unrecognized filetype of ''%s''', filename));
end

% set up the axes of the volume in voxel coordinates
nx = size(img,1);
ny = size(img,2);
nz = size(img,3);
mri.dim = [nx ny nz];
% store the anatomical data
mri.anatomy = img;
% store the header with all fileformat specific details
mri.hdr = hdr;
try
  % if present, store the homogenous transformation matrix
  mri.transform = transform;
end
try
  % try to add units
  mri = ft_convert_units(mri);
end

author: PatternRecognition

repo: OpenBMI


text:
function [ y, j ] = j_borrow_gregorian ( y, j )

%*****************************************************************************80
%
%% J_BORROW_GREGORIAN borrows year-days from years in a Gregorian date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, a YJ date.
%
  while ( j <= 0 )

    y = y - 1;

    days = year_length_gregorian ( y );

    j = j + days;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = f_02 ( x )

%*****************************************************************************80
%
%% F_02 evaluates 2*x-exp(-x).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the point at which F is to be evaluated.
%
%    Output, real VALUE, the value of the function at X.
%
  value = 2.0 * x - exp ( - x );

  return
end


text:
function p = detect_quadratic_disjoints(p)

p.quadraticdisjoints = [];
if any(p.K.l)
    top = startofLPCone(p.K);
    for i = 1:p.K.l
        row = p.F_struc(top,:);
        if nnz(row)<=3
            a = row(1);
            [~,var,val] = find(row(2:end));
            if length(var)==1
                if p.variabletype(var)==2
                    if val > 0 && a < 0
                        x = find(p.monomtable(var,:));
                        p.quadraticdisjoints = [p.quadraticdisjoints [x;sqrt(-a);-sqrt(-a)]];
                        % val*x^2 - a >= 0
                        % i.e >=sqrt(a) or <=-sqrt(a)
                    end
                end
            end
        end
        top = top + 1;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


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

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

% Last Modified by GUIDE v2.5 09-Sep-2013 22:52:08

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

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

% --- Executes during object creation, after setting all properties.
function f_CreateFcn(hObject, eventdata, handles)
% hObject    handle to f (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function f_Callback(hObject, eventdata, handles)
% hObject    handle to f (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of f as text
%        str2double(get(hObject,'String')) returns contents of f as a double

% --- Executes on button press in plot.
function plot_Callback(hObject, eventdata, handles)
syms x;
fun=(x.^2)-1;
ezplot(fun);
% hObject    handle to plot (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes during object creation, after setting all properties.
function x0_CreateFcn(hObject, eventdata, handles)
% hObject    handle to x0 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function x0_Callback(hObject, eventdata, handles)
% hObject    handle to x0 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of x0 as text
%        str2double(get(hObject,'String')) returns contents of x0 as a double

% --- Executes during object creation, after setting all properties.
function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double

% --- Executes during object creation, after setting all properties.
function edit5_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function edit5_Callback(hObject, eventdata, handles)
% hObject    handle to edit5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit5 as text
%        str2double(get(hObject,'String')) returns contents of edit5 as a double

% --- Executes on button press in solve.
function solve_Callback(hObject, eventdata, handles)

xi=str2num(get(handles.x0,'string'));

for i=1:20

xn=xi-(polyval([1 0 -1],[xi])/polyval(polyder([1 0 -1]),[xi]));
disp(xn);
set(handles.r,'string',xn);
xi=xn;
i=i+1;

end
% hObject    handle to solve (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

author: Sable

repo: mcbench-benchmarks


text:
colors =  ['r' 'b' 'k' 'g' 'c' 'y' 'm' ...
	   'r' 'b' 'k' 'g' 'c' 'y' 'm'];
symbols = ['o' 'x' '+' '>' '<' '^' 'v' ...
	   '*' 'p' 'h' 's' 'd' 'o' 'x'];
for i=1:length(colors)
  styles{i} = sprintf('-%s%s', colors(i), symbols(i));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [nVar, lowerbound, upperbound] = getVarInfo(obj)
    % The function returns the dimension, upper/lower limits of NLP
    % variables.
    %
    % Return values:
    %  nVar: the total dimension of all NLP variables @type
    %  integer
    %  lowerbound: the lower limits @type colvec
    %  upperbound: the upper limits @type colvec

    assert(~isempty(obj.VariableArray),...
        'NonlinearProgram:emptyVarArray',...
        ['No variable has been definied.\n',...
        'Please define NLP variables first.\n']);

    nVar = sum([obj.VariableArray.Dimension]);
    lowerbound = vertcat(obj.VariableArray.LowerBound);
    upperbound = vertcat(obj.VariableArray.UpperBound);

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function initRegularisation(solver,varargin)

epsilon = get_option(varargin,'epsilon',solver.psi.halfwidth*2);

% a discrete Laplacian
A = -dot_outer(solver.S3G,solver.S3G,'epsilon',epsilon);
d = full(sum(A)).';

solver.RM = solver.lambda * spdiags(-1-d,0,A);

author: mtex-toolbox

repo: mtex

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


text:
function window_sz = get_search_window_test( target_sz, im_sz)
% GET_SEARCH_WINDOW

% if(target_sz(1)/target_sz(2) > 2)
%     % For objects with large height, we restrict the search window with padding.height
%     window_sz = floor(target_sz.*[1+padding.height, 1+padding.generic]);
%
% elseif(prod(target_sz)/prod(im_sz(1:2)) > 0.05)
%     % For objects with large height and width and accounting for at least 10 percent of the whole image,
%     % we only search 2x height and width
%     window_sz=floor(target_sz*(1+padding.large));
%
% else
%     %otherwise, we use the padding configuration
%     window_sz = floor(target_sz * (1 + padding.generic));

ratio=target_sz(1)/target_sz(2);
if ratio>1
    window_sz=round(target_sz.*[2,2*ratio]);
else
    window_sz=round(target_sz.*[2/ratio,2]);
end
%
% %window_sz=round(target_sz.*[5,9]);
window_sz=window_sz-mod(window_sz,2)+1;

end

repo: CREST-Release


text:
function [ a, irow, icol ] = i4mat_red ( m, n, a )

%*****************************************************************************80
%
%% I4MAT_RED divides out common factors in a row or column of an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in the matrix.
%
%    Input, integer N, the number of columns in the matrix.
%
%    Input, integer A(M,N), the matrix to be reduced.
%
%    Output, integer A(M,N), the reduced matrix.  The greatest common factor
%    in any row or column is 1.
%
%    Output, integer IROW(M), the row factors that were divided out.
%
%    Output, integer ICOL(N), the column factors that were divided out.
%
  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_RED - Warning!\n' );
    fprintf ( 1, '  M must be greater than 0.\n' );
    fprintf ( 1, '  Input M = %d\n', m );
    error ( 'I4MAT_RED - Fatal error!' );
  end

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_RED - Warning!\n' );
    fprintf ( 1, '  N must be greater than 0.\n' );
    fprintf ( 1, '  Input N = %d\n', n );
    error ( 'I4MAT_RED - Fatal error!' );
  end
%
%  Remove factors common to a column.
%
  incx = 1;
  for j = 1 : n
    [ a(1:m,j), icol(j) ] = i4vec_red ( m, a(1:m,j), incx );
  end
%
%  Remove factors common to a row.
%
  incx = 1;
  for i = 1 : m
    [ a(i,1:n), irow(i) ] = i4vec_red ( n, a(i,1:n), incx );
  end

  return
end


text:
function cdf = gumbel_cdf ( x )

%*****************************************************************************80
%
%% GUMBEL_CDF evaluates the Gumbel CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Output, real CDF, the value of the CDF.
%
  cdf = exp ( - exp ( - x ) );

  return
end


text:
function test6
%TEST6 test for BTF
% Requires UFget
% Example:
%   test6
% See also btf, maxtrans, strongcomp, dmperm, UFget,
%   test1, test2, test3, test4, test5.

% Copyright 2007, Timothy A. Davis, University of Florida

quick2 = [ ...
 1522 -272  1463  1521   460 1507  -838 1533 -1533 -1456 -1512   734   211 ...
 -385 -735   394  -397  1109 -744  ...
 -734 -375 -1200 -1536  -837  519  -519  520  -520   189  -189   454   385 ...
  387 -387   384  -384   386 -386   388 -388   525  -525   526  -526   735 ...
 1508  209   210  1243 -1243 1534  -840 1234 -1234   390  -390   392  -392 ...
 -394 1472  1242 -1242   389 -389   391 -391   393  -393  1215 -1215  1216 ...
-1216  736  -736   737  -737  455  -455 -224  -839  1426 -1426 -1473   396 ...
 -396  398  -398   400  -400  402  -402  404  -404 -1531   395  -395   397 ...
  399 -399   401  -401   403 -403   405 -405  -738  -739  1459 -1459  1111 ...
 1110  376  -376   284  -284 -740  -742 -741  -743  1293 -1293   452   920 ...
 -745 -446  1462 -1461   448 -448   283 -283  1502 -1502  1292 -1292  1503 ...
-1503 1291 -1291   445  -445 -746  -747 1300 -1300   435  -435 -1343 -1345 ...
-1344 1305 -1305   921 -1513 1307 -1307 1369 -1369  1374 -1374  1377 ...
-1377  748  -748  -749  1510  922  -922 ] ;

index = UFget ;
nmat = length (quick2) ;
dopause = 0 ;

h = waitbar (0, 'BTF test 6 of 6') ;

try

    for k = 1:nmat

        waitbar (k/nmat, h) ;

        i = quick2 (k) ;
        Prob = UFget (abs (i), index) ;
        disp (Prob) ;
        if (i < 0)
            fprintf ('transposed\n') ;
            A = Prob.A' ;
            [m n] = size (A) ;
            if (m == n)
                if (nnz (spones (A) - spones (Prob.A)) == 0)
                    fprintf ('skip...\n') ;
                    continue ;
                end
            end
        else
            A = Prob.A ;
        end

        tic
        [p1,q1,r1,work1] = btf (A) ;
        t1 = toc ;
        n1 = length (r1) - 1 ;
        m1 = nnz (diag (A (p1, abs (q1)))) ;

        limit = work1/nnz(A) ;

        fprintf ('full search: %g * nnz(A)\n', limit) ;

        works = linspace(0,limit,9) ;
        works (1) = eps ;
        nw = length (works) ;

        T2 = zeros (nw, 1) ;
        N2 = zeros (nw, 1) ;
        M2 = zeros (nw, 1) ;

        T2 (end) = t1 ;
        N2 (end) = n1 ;
        M2 (end) = m1 ;

        fprintf ('full time %10.4f   blocks %8d  nnz(diag) %8d\n\n', t1, n1, m1) ;

        subplot (3,4,4) ;
        drawbtf (A, p1, abs (q1), r1) ;
        title (Prob.name, 'Interpreter', 'none') ;

        for j = 1:nw-1

            maxwork = works (j) ;

            tic
            [p2,q2,r2,work2] = btf (A, maxwork) ;
            t2 = toc ;
            n2 = length (r2) - 1 ;
            m2 = nnz (diag (A (p2, abs (q2)))) ;
            T2 (j) = t2 ;
            N2 (j) = n2 ;
            M2 (j) = m2 ;

            fprintf ('%9.1f %10.4f   blocks %8d  nnz(diag) %8d\n', ...
                maxwork, t2, n2, m2) ;

            subplot (3,4,4+j) ;
            drawbtf (A, p2, abs (q2), r2) ;
            title (sprintf ('%g', maxwork)) ;

            ss = [1:j nw] ;

            subplot (3,4,1) ;
            plot (works(ss), T2(ss), 'o-') ;  title ('time vs work') ;
            axis ([0 limit 0 max(0.1,max(T2))]) ;

            subplot (3,4,2) ;
            plot (works(ss), N2(ss), 'o-') ; title ('blocks vs work') ;
            axis ([0 limit 0 n1]) ;

            subplot (3,4,3) ;
            plot (works(ss), M2(ss), 'o-') ; title ('nnz(diag) vs work') ;
            axis ([0 limit 0 m1]) ;
            drawnow

        end
        fprintf ('full time %10.4f   blocks %8d  nnz(diag) %8d\n', t1, n1, m1) ;

        if (dopause)
            input ('hit enter: ') ;
        end

    end

catch
    % out-of-memory is OK, other errors are not
    disp (lasterr) ;
    if (isempty (strfind (lasterr, 'Out of memory')))
        error (lasterr) ;                                                   %#ok
    else
        fprintf ('test terminated early, but otherwise OK\n') ;
    end
end

close (h) ;


text:
function [s_true,e_spat,e_interf,e_artif]=bss_decomp_mtifilt(se,S,j,L)

% BSS_DECOMP_MTIFILT Decomposition of an estimated source image into four
% components representing respectively the true source image, spatial (or
% filtering) distortion, interference and artifacts, derived from the true
% source images using multichannel time-invariant filters.
%
% [s_true,e_spat,e_interf,e_artif]=bss_decomp_mtifilt(se,S,j,L)
%
% Inputs:
% se: I x T matrix containing the estimated source image (one row per channel)
% S: I x T x J matrix containing the true source images
% j: source index corresponding to the estimated source image in S
% L: length of the multichannel time-invariant filters in samples
%
% Outputs:
% s_true: I x T matrix containing the true source image, i.e. S(:,:,j) (one row per channel)
% e_spat: I x T matrix containing the spatial (or filtering) distortion component
% e_interf: I x T matrix containing the interference component
% e_artif: I x T matrix containing the artifacts component

%%% Errors %%%
if nargin<4, error('Not enough input arguments.'); end
[Ie,Te]=size(se);
[I,T,J]=size(S);
if I~=Ie, error('The number of channels of the true source images and the estimated source image must be equal.'); end
if T~=Te, error('The duration of the true source images and the estimated source image must be equal.'); end

%%% Decomposition %%%
% True source image
s_true=[S(:,:,j),zeros(I,L-1)];
% Spatial (or filtering) distortion
e_spat=project(se,S(:,:,j),L)-s_true;
% Interference
e_interf=project(se,S,L)-s_true-e_spat;
% Artifacts
e_artif=[se,zeros(I,L-1)]-s_true-e_spat-e_interf;

return;

function sproj=project(se,S,L)

% Least-squares projection of each channel of se on the subspace spanned by
% delayed versions of the channels of S, with delays between 0 and L-1

[I,T,J]=size(S);
S=reshape(permute(S,[1 3 2]),I*J,T);

%%% Computing coefficients of least squares problem via FFT %%%
% Zero padding and FFT of input data
S=[S,zeros(I*J,L-1)];
se=[se,zeros(I,L-1)];
N=2^nextpow2(T+L-1);
Sf=fft(S,N,2);
sef=fft(se,N,2);
% Inner products between delayed versions of S
G=zeros(I*J*L);
for k1=0:I*J-1
    for k2=0:k1
        SSf=Sf(k1+1,:).*conj(Sf(k2+1,:));
        SSf=real(ifft(SSf));
        SS=toeplitz(SSf([1 N:-1:N-L+2]),SSf(1:L));
        G(k1*L+1:k1*L+L,k2*L+1:k2*L+L)=SS;
        G(k2*L+1:k2*L+L,k1*L+1:k1*L+L)=SS.';
    end
end
% Inner products between se and delayed versions of S
D=zeros(I*J*L,I);
for k=0:I*J-1
    for i=1:I
        Ssef=Sf(k+1,:).*conj(sef(i,:));
        Ssef=real(ifft(Ssef,[],2));
        D(k*L+1:k*L+L,i)=Ssef(:,[1 N:-1:N-L+2]).';
    end
end

%%% Computing projection %%%
% Distortion filters
C=G\D;
C=reshape(C,L,I*J,I);
% Filtering
sproj=zeros(I,T+L-1);
for k=1:I*J
    for i=1:I
        sproj(i,:)=sproj(i,:)+fftfilt(C(:,k,i).',S(k,:));
    end
end

return;

repo: Voice-Separation-and-Enhancement


text:
function [str, x_bytes, unit] = bytes2str(bytes, precision)
%BYTES2STR Returns the number of bytes in a more readable format.
% Usage:
%   bytes2str(bytes)
%   str = bytes2str(bytes)
%   str = bytes2str(filename)
%   str = bytes2str(_, precision)
%   [str, x_bytes, unit] = bytes2str(...)

if ischar(bytes)
    bytes = get_filesize(bytes);
end
if nargin < 2
    precision = 3;
end

units = {'bytes', 'KB', 'MB', 'GB', 'TB', 'PB'};

% Find closest units
base = floor(log(bytes) / log(1024));

% Convert to new units
x_bytes = bytes * (1024 ^ -base);

% Convert to string
unit = units{base + 1};
str = sprintf('%s %s', num2str(x_bytes, precision), unit);

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function handleKey(dummy1,dummy2)
global finish
finish = 1;
fprintf('Execution interrupted by keypress.\n');
end % by pressing any key, the process will exit

author: flyers

repo: drone-tracking


text:
function [images images2 images3 view]= RawToCommonAtlas(images, corners, safirData, safirResult, view, flip_lr)
% From the results of atlas fitting, get the data morphed into an individual
% atlas space (images2), and that morphed into a common atlas space (images3)
%
% [images images2 images3 view]= RawToCommonAtlas(images, corners, safirData, safirResult, [view], [flip_lr])
%
% INPUT
%  images: a structure including raw data (M1,M2,M3) and inital atlas (A1,A2)
%  corners: corners of the initial atlas
%  safirData: data structure for the atlas fitting GUI
%  safirResult: morphing matrix
%  view: has to be specified if you want to get ROIs
%  flip_lr: flip
%
% OUTPUT
%  images: raw data with the fitted atlas
%  images2: the data morphed into an individual atlas space
%  images3: the data morphed into a common atlas space
%  view: has to be specified if you want to get ROIs
% Example
%  [images images2 images3]= RawToCommonAtlas(images, corners, 'safirData', 'safirResult');
%
% 07/03 KA wrote it

% Check arguments
if notDefined('images'), error('images should be defined'); end
if notDefined('corners'), error('corners should be defined'); end
if notDefined('safirData'), safirData = 'safirData'; end
if notDefined('safirResult'), safirData = 'safirResult'; end
if notDefined('flip_lr'), flip_lr=0; end

load(safirResult)
load(safirData,'maxLevel')

if size(corners,2)==6 % LO,TO
    area_num=4;
elseif size(corners,2)==5 % V1-3
    area_num=5;
end

yOpt_inv = inverse_yOpt(yOpt,maxLevel);
if isfield(images,'A3')
    [dM1m dM2m dM3m Edges dEdges A1m A2m M1m M2m M3m W1m W2m COm dW1m dW2m dCOm dAreasImgm dA1m dA2m dA3m] = showResults2(yOpt,yOpt_inv,corners, safirData,images);
else
    [dM1m dM2m dM3m Edges dEdges A1m A2m M1m M2m M3m W1m W2m COm dW1m dW2m dCOm dAreasImgm dA1m dA2m] = showResults2(yOpt,yOpt_inv,corners, safirData,images);
end

images_size = size(images.A1,1);
if strfind(safirData,'left')>0
    curSlice = 1;
elseif  strfind(safirData,'right')>0
    curSlice = 2;
end
if notDefined('view')==0
    view = transformAtlasToROIs(round(dAreasImgm), images_size, curSlice, view);
end

% raw data spcae
images.M1 = M1m;
images.M2 = M2m;
images.M3 = M3m;
images.W1 = W1m;
images.W2 = W2m;
images.CO = COm;
images.dA1 = dA1m;
images.dA2 = dA2m;
if isfield(images,'A3')
    images.dA3 = dA3m;
end
images.dEdges = dEdges;
images.dareasImg = round(dAreasImgm);

% individual atlas space
images2.A1 = A1m;
images2.A2 = A2m;
images2.M1 = dM1m;
images2.M2 = dM2m;
images2.M3 = dM3m;
images2.W1 = dW1m;
images2.W2 = dW2m;
images2.CO = dCOm;
images2.areasImg = round(dAreasImgm);

% For flipped image (since sometimes flat view itself is flipped)
if flip_lr
    images2.M1 = fliplr(images2.M1);
    images2.M2 = fliplr(images2.M2);
    images2.M3 = fliplr(images2.M3);
    images2.A1 = fliplr(images2.A1);
    images2.A2 = fliplr(images2.A2);
    images2.W1 = fliplr(images2.W1);
    images2.W2 = fliplr(images2.W2);
    images2.CO = fliplr(images2.CO);
    images2.areasImg = fliplr(images2.areasImg);

    for i=1:size(corners,2)
        corners{i}(:,1) = size(images.A1,2)+1-corners{i}(:,1);
    end
end

% scale corner coordinates
for i=1:size(corners,2)
    corners_scaled{i} = corners{i}*(size(images.A1,1)-1)/(size(images.A1,1)-1);
%     corners_scaled{i} = corners{i};
end

if size(corners,2)==6 % LO-1/2, TO-1/2
    images3.A1 = transform_image_LO_TO(images2.A1,corners,'bicubic');
    images3.A2 = transform_image_LO_TO(images2.A2,corners,'bicubic');
    images3.M1 = transform_image_LO_TO(images2.M1,corners,'bicubic');
    images3.M2 = transform_image_LO_TO(images2.M2,corners,'bicubic');
    images3.M3 = transform_image_LO_TO(images2.M3,corners,'bicubic');
    images3.W1 = transform_image_LO_TO(images2.W1,corners,'bicubic');
    images3.W2 = transform_image_LO_TO(images2.W2,corners,'bicubic');
    images3.CO = transform_image_LO_TO(images2.CO,corners,'bicubic');
    images3.CO2=images3.CO;
    images3.CO2(find(images3.CO2<0.1))=0;
    images3.COmask = double(images3.CO>=0.1);
    images3.areasImg = round(transform_image_LO_TO(images2.areasImg,corners,'nearest'));
elseif size(corners,2)==5 % V1-3
    images3.A1 = transform_image_V1_V3(images2.A1,corners,'bicubic');
    images3.A2 = transform_image_V1_V3(images2.A2,corners,'bicubic');
    images3.M1 = transform_image_V1_V3(images2.M1,corners,'bicubic');
    images3.M2 = transform_image_V1_V3(images2.M2,corners,'bicubic');
    images3.M3 = transform_image_V1_V3(images2.M3,corners,'bicubic');
    images3.W1 = transform_image_V1_V3(images2.W1,corners,'bicubic');
    images3.W2 = transform_image_V1_V3(images2.W2,corners,'bicubic');
    images3.CO = transform_image_V1_V3(images2.CO,corners,'bicubic');
    images3.CO2=images3.CO;
    images3.CO2(find(images3.CO2<0.1))=0;
    images3.COmask = double(images3.CO>=0.1);
    images3.areasImg = round(transform_image_V1_V3(images2.areasImg,corners,'nearest'));
end

if ~notDefined('images3')
    images3.M1(find(isnan(images3.CO)))=0;
    images3.M2(find(isnan(images3.CO)))=0;
    images3.M3(find(isnan(images3.CO)))=0;
    images3.CO(find(isnan(images3.CO)))=0;

    images3.M1(find(images3.CO>1))=0;
    images3.M2(find(images3.CO>1))=0;
    images3.M3(find(images3.CO>1))=0;
    images3.CO(find(images3.CO>1))=0;
end

if size(corners,2)==6 % LO-1/2, TO-1/2
%     corners2 =  [30 40; 90 40; 90 80; 30 80; 30 40; 45 40; 45 80; 60 80; 60 40; 75 40; 75 80];
    corners2 =  [30 40; 105 40; 105 80; 105 40; 90 40; 90 80; 30 80; 30 40; 45 40; 45 80; 60 80; 60 40; 75 40; 75 80; 90 80; 105 80];
else % V1-3
    corners2 =  [15 40; 105 40; 105 80; 15 80; 15 40; 30 40; 30 80; 45 80; 45 40; 75 40; 75 80; 90 80; 90 40];
end

tmp=load('WedgeMapLeft');
figure
subplot(331)
imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(images.M2.*double(images.CO>0.1),dEdges,[tmp.modeInformation.cmap(119:128,:);tmp.modeInformation.cmap(65:118,:)],2*pi))
% contourf(images.M2.*double(images.CO>0.1))
set(gca,'YDir','reverse');
axis equal
axis([1 images_size 1 images_size])
title('Angle map with fitted atlas')
set(gca,'XTickLabel','')
set(gca,'YTickLabel','')
axis off

subplot(332)
hold on
imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(images2.M2.*double(images2.CO>0.1),[],[tmp.modeInformation.cmap(119:128,:);tmp.modeInformation.cmap(65:118,:)],2*pi))
for ii=1:size(corners,2)
    plot([corners{ii}(:,1);corners{ii}(1,1)],[corners{ii}(:,2);corners{ii}(1,2)],'k','LineWidth',2)
end
if size(corners,2)==5
    plot([corners{5}(:,1);corners{5}(1,1)],[corners{5}(:,2);corners{5}(1,2)],'k','LineWidth',2)
end
% plot([corners{7}(:,1);corners{7}(1,1)],[corners{7}(:,2);corners{7}(1,2)],'k','LineWidth',2)
set(gca,'YDir','reverse');
axis equal
axis([1 images_size 1 images_size])
title('Angle (raw data space)')
title('Deformed angle map with initial atlas')
set(gca,'XTickLabel','')
set(gca,'YTickLabel','')
axis off

if ~notDefined('images3')
    subplot(333)
    hold on
    imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(images3.M2.*double(images3.CO>0.1),[],[tmp.modeInformation.cmap(119:128,:);tmp.modeInformation.cmap(65:118,:)],2*pi))
    plot(corners2(:,1),corners2(:,2),'k','LineWidth',2)
    set(gca,'YDir','reverse');
    axis equal
    axis([1 120 1 120])
    title('Deformed angle map with common atlas')
    set(gca,'XTickLabel','')
    set(gca,'YTickLabel','')
    axis off
end

subplot(334)
imagesc(0:images_size/127:images_size,0:images_size/127:images_size,mergedImage(images.M1.*double(images.CO>0.1),dEdges,hsvTbCmap(0,256), 15))
set(gca,'YDir','reverse');
axis equal
axis([1 images_size 1 images_size])
title('Eccentricity map with fitted atlas')
set(gca,'XTickLabel','')
set(gca,'YTickLabel','')
axis off

subplot(335)
hold on
imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(images2.M1.*double(images2.CO>0.1),[],hsvTbCmap(0,256),15))
for ii=1:size(corners,2)
    plot([corners{ii}(:,1);corners{ii}(1,1)],[corners{ii}(:,2);corners{ii}(1,2)],'k','LineWidth',2)
end
set(gca,'YDir','reverse');
axis equal
axis([1 images_size 1 images_size])
title('Deformed eccentricity map with initial atlas')
set(gca,'XTickLabel','')
set(gca,'YTickLabel','')
axis off

if ~notDefined('images3')
    subplot(336)
    hold on
    imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(images3.M1.*double(images3.CO>0.1),[],hsvTbCmap(0,256),15))
    plot(corners2(:,1),corners2(:,2),'k','LineWidth',2)
    set(gca,'YDir','reverse');
    axis equal
    axis([1 120 1 120])
    title('Deformed eccentricity map with common atlas')
    set(gca,'XTickLabel','')
    set(gca,'YTickLabel','')
    axis off
end

tmp=cool_springCmap;
subplot(337)
imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(images.M3.*double(images.CO>0.1),dEdges,tmp(129:224,:), 15))
set(gca,'YDir','reverse');
axis equal
axis([1 images_size 1 images_size])
title('pRF size map with fitted atlas')
set(gca,'XTickLabel','')
set(gca,'YTickLabel','')
axis off

subplot(338)
hold on
imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(images2.M3.*double(images2.CO>0.1),[],tmp(129:224,:),15))
for ii=1:size(corners,2)
    plot([corners{ii}(:,1);corners{ii}(1,1)],[corners{ii}(:,2);corners{ii}(1,2)],'k','LineWidth',2)
end
set(gca,'YDir','reverse');
axis equal
axis([1 images_size 1 images_size])
title('Deformed pRF size map with initial atlas')
set(gca,'XTickLabel','')
set(gca,'YTickLabel','')
axis off

if ~notDefined('images3')
    subplot(339)
    hold on
    imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(images3.M3.*double(images3.CO>0.1),[],tmp(129:224,:),15))
    plot(corners2(:,1),corners2(:,2),'k','LineWidth',2)
    set(gca,'YDir','reverse');
    axis equal
    axis([1 120 1 120])
    title('Deformed pRF size map with common atlas')
    set(gca,'XTickLabel','')
    set(gca,'YTickLabel','')
    axis off

    % mask = (images3.CO>0.1).*(images3.M2>pi/2-pi/4).*(images3.M2<pi/2*3+pi/4);
    mask = double(images3.CO>0.1).*double(double(images3.CO<=1));
    mask_nan = mask;
    mask_nan(find(mask==0))=nan;

    % Angle modulation averaged across iso-ecc lines
    polar_ave=nan(3,size(round(10/(images_size-1)*(2^maxLevel-1)):round(110/(images_size-1)*(2^maxLevel-1)),2));
    polar_min=nan(3,size(round(10/(images_size-1)*(2^maxLevel-1)):round(110/(images_size-1)*(2^maxLevel-1)),2));
    polar_max=nan(3,size(round(10/(images_size-1)*(2^maxLevel-1)):round(110/(images_size-1)*(2^maxLevel-1)),2));
    for i=round(10/(images_size-1)*(2^maxLevel-1)):round(110/(images_size-1)*(2^maxLevel-1))
        for k=1:3
            tmp=[];
            for j=round((80-k*10)/(images_size-1)*(2^maxLevel-1)):round((90-k*10)/(images_size-1)*(2^maxLevel-1))
                % select the data for the averaging
                if mask(j,i)==1
                    tmp = [tmp [images3.M2(j,i);images3.CO(j,i)]];
                end
            end
            if sum(size(tmp))>0
                polar_ave(k,i-round(10/(images_size-1)*(2^maxLevel-1))+1)=sum(tmp(1,:).*tmp(2,:))/sum(tmp(2,:));
                polar_min(k,i-round(10/(images_size-1)*(2^maxLevel-1))+1)=min(tmp(1,:));
                polar_max(k,i-round(10/(images_size-1)*(2^maxLevel-1))+1)=max(tmp(1,:));
            end
        end
    end

    % Ecc modulation averaged across iso-polar lines
    ecc_ave=nan(area_num,44);
    for j=round(40/(images_size-1)*(2^maxLevel-1)):round(80/(images_size-1)*(2^maxLevel-1))
        for k=1:area_num
            tmp=[];
            for i=round((15+k*15)/(images_size-1)*(2^maxLevel-1)):round((30+k*15)/(images_size-1)*(2^maxLevel-1))
                if mask(j,i)==1
                    tmp = [tmp [images3.M1(j,i);images3.CO(j,i)]];
                end
            end
            if sum(size(tmp))>0
                ecc_ave(k,round(80/(images_size-1)*(2^maxLevel-1))+1-j)=sum(tmp(1,:).*tmp(2,:))/sum(tmp(2,:));
            end
        end
    end

    upper_model = [pi:pi/2/16:pi*1.5 pi*1.5-pi/2/16:-pi/2/16:pi pi+pi/2/16:pi/2/16:pi*1.5 pi*1.5-pi/2/16:-pi/2/16:pi];
    lower_model = [pi:pi/2/16:pi*1.5 pi*1.5-pi/2/16:-pi/2/16:pi pi+pi/2/16:pi/2/16:pi*1.5 pi*1.5-pi/2/16:-pi/2/16:pi]-pi*0.5;
    hemi_model = [pi*0.5:pi/16:pi*1.5 pi*1.5-pi/16:-pi/16:pi*0.5 pi*0.5+pi/16:pi/16:pi*1.5 pi*1.5-pi/16:-pi/16:pi*0.5];

    % Phase plot along iso-ecc lines and ecc plot along iso-polar lines
    figure
    for k=1:3
        subplot(2,3,k)
        hold on
        plot(-20:100/(round(110/(images_size-1)*(2^maxLevel-1))-round(10/(images_size-1)*(2^maxLevel-1))):80, ...
            images3.M2(round((80-k*10)/(images_size-1)*(2^maxLevel-1)):round((90-k*10)/(images_size-1)*(2^maxLevel-1)),round(10/(images_size-1)*(2^maxLevel-1)):round(110/(images_size-1)*(2^maxLevel-1)))' ...
            .*mask_nan(round((80-k*10)/(images_size-1)*(2^maxLevel-1)):round((90-k*10)/(images_size-1)*(2^maxLevel-1)),round(10/(images_size-1)*(2^maxLevel-1)):round(110/(images_size-1)*(2^maxLevel-1)))'-pi,'k','LineWidth',1)
        plot(-20:100/(round(110/(images_size-1)*(2^maxLevel-1))-round(10/(images_size-1)*(2^maxLevel-1))):80, polar_ave(k,:)-pi,'r','LineWidth',2)
        %     plot(0:70/74:70, polar_min(k,:)-pi,'r:','LineWidth',2)
        %     plot(0:70/74:70, polar_max(k,:)-pi,'r:','LineWidth',2)
        %     plot(0:80/85:80, upper_model-pi,'b','LineWidth',2);
        %     plot(0:80/85:80, lower_model-pi,'g','LineWidth',2);
        %     plot(0:80/85:80, hemi_model-pi,'m','LineWidth',2);
        xlabel('Atlas distance')
        ylabel('Angle from HM')
        if size(corners,2)==6 % LO,TO
            xlim([-10 70])
        else
            xlim([-20 80])
        end
        ylim([-2 2])
        if size(corners,2)==6 % LO,TO
            text(3,1.5,'LO1')
            text(18,1.5,'LO2')
            text(33,1.5,'TO1')
            text(48,1.5,'TO2')
        else
            text(-10,1.8,'V3v')
            text(5,1.8,'V2v')
            text(28,1.8,'V1')
            text(50,1.8,'V2d')
            text(65,1.8,'V3d')
        end

        if k==1
            title('Ecc = 0-4 deg')
        elseif k==2
            title('Ecc = 4-8 deg')
        else
            title('Ecc = 8-12 deg')
        end
    end

    for k=1:area_num
        subplot(2,area_num,area_num+k)
        hold on
        plot(40:-40/(round(80/(images_size-1)*(2^maxLevel-1))-round(40/(images_size-1)*(2^maxLevel-1))):0, ...
            images3.M1(round(40/(images_size-1)*(2^maxLevel-1)):round(80/(images_size-1)*(2^maxLevel-1)),round((15+15*k)/(images_size-1)*(2^maxLevel-1)):round((30+15*k)/(images_size-1)*(2^maxLevel-1))) ...
            .*mask_nan(round(40/(images_size-1)*(2^maxLevel-1)):round(80/(images_size-1)*(2^maxLevel-1)),round((15+15*k)/(images_size-1)*(2^maxLevel-1)):round((30+15*k)/(images_size-1)*(2^maxLevel-1))),'k','LineWidth',1)
        plot(0:(round(80/(images_size-1)*(2^maxLevel-1))-round(40/(images_size-1)*(2^maxLevel-1))):40, ecc_ave(k,:),'r','LineWidth',2)
        xlabel('Atlas distance')
        ylabel('Eccentricity (deg)')
        xlim([0 40])
        ylim([0 15])
        if size(corners,2)==6 % LO,TO
            if k==1
                title('LO1')
            elseif k==2
                title('LO2')
            elseif k==3
                title('TO1')
            else
                title('TO2')
            end
        elseif size(corners,2)==5 % V1-3
            if k==1
                title('V3v')
            elseif k==2
                title('V2v')
            elseif k==3
                title('V1')
            elseif k==4
                title('V2d')
            else
                title('V3d')
            end
        end
    end
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [modelSize] = modelSizes(model)
%%%
%
% This method prints out detailed sizes of individual matrices and return
% the total model size.
% model is a struct, in which each element could be a matrix or a cell of
% matrices.
%
% Thang Luong @ 2015, <lmthang@stanford.edu>
%
%%%
  fieldNames = fields(model);

  % compute model size
  modelSize = 0;
  for ii=1:length(fieldNames)
    field = fieldNames{ii};
    if iscell(model.(field))
      for jj=1:length(model.(field))
        modelSize = modelSize + numel(model.(field){jj});
        modelSizes.(field){jj} = numel(model.(field){jj});
      end
    else
      modelSize = modelSize + numel(model.(field));
      modelSizes.(field) = numel(model.(field));
    end
  end
  fprintf(2, '  Model size = %d, individual sizes: %s\n', modelSize, wInfo(modelSizes, 0));
end

author: lmthang


text:
function [m, symbol, filterind] = model_add_terminal(m, varargin)
% Add a filter to the model.
%   [m, symbol, filterind] = model_add_terminal(m, varargin)
%
% Return values
%   m         Updated model
%   symbol    Newly created terminal symbol
%   filterind Index of the newly created filter associated with the terminal
%
% Arguments
%   m         Model to update
%   varargin  (key, value) pairs that can specify the following:
%   key                         value
%   ---                         -----
%   w                           Filter coefficients
%   flip                        True or false (default)
%   blocklabel                  model.blocks index
%   mirror_terminal             Terminal symbol to horizontally mirror

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

valid_opts = {'w', 'blocklabel', 'flip', 'mirror_terminal'};
opts = getopts(varargin, valid_opts);

if opts.isKey('mirror_terminal')
  src_terminal = opts('mirror_terminal');
  fi = m.symbols(src_terminal).filter;
  opts('blocklabel') = m.filters(fi).blocklabel;
  opts('w')          = flipfeat(model_get_block(m, m.filters(fi)));
  opts('flip')       = ~m.filters(fi).flip;
end

if opts.isKey('w')
  w = opts('w');
else
  error('argument ''w'' is required');
end

if opts.isKey('blocklabel')
  blocklabel = opts('blocklabel');
else
  [m, blocklabel] = model_add_block(m, ...
                                    'type', block_types.Filter, ...
                                    'w', w);
end

if opts.isKey('flip')
  flip = opts('flip');
else
  flip = false;
end

% get index for new filter
j = m.numfilters + 1;
m.numfilters = j;

m.filters(j).blocklabel = blocklabel;
m.filters(j).size       = [size(w, 1) size(w, 2)];
m.filters(j).flip       = flip;

% new symbol for terminal associated with filter f
[m, i] = model_add_symbol(m, 'T');
m.symbols(i).filter = j;
m.filters(j).symbol = i;

filterind = j;
symbol = i;

author: rbgirshick

repo: voc-dpm

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


text:
function [tx, ty] = orientedBoxToPolygon(obox)
%ORIENTEDBOXTOPOLYGON Convert an oriented box to a polygon (set of vertices)
%
%   POLY = orientedBoxToPolygon(OBOX);
%   Converts the oriented box OBOX given either as [XC YC W H] or as
%   [XC YC W H THETA] into a 4-by-2 array of double, containing coordinates
%   of box vertices.
%   XC and YC are center of the box. W and H are the width and the height
%   (dimension in the main directions), and THETA is the orientation, in
%   degrees between 0 and 360.
%
%   Example
%     OBOX = [20 10  40 20 0];
%     RECT = orientedBoxToPolygon(OBOX)
%     RECT =
%         -20 -10
%          20 -10
%          20  10
%         -20  10
%
%
%   See also:
%   polygons2d, orientedBox, drawOrientedBox, rectToPolygon
%

%   ---------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% INRA - TPV URPOI - BIA IMASTE
% created the 06/04/2005.
%

%   HISTORY
%   2011-10-09 rewrite from rectAsPolygon to orientedBoxToPolygon
%   2016: Simplify by JuanPi Carbajal

% extract box parameters
theta = 0;
x = obox(1);
y = obox(2);
w = obox(3) / 2;  % easier to compute with w and h divided by 2
h = obox(4) / 2;
if length(obox) > 4
    theta = obox(5);
end

v = [cosd(theta); sind(theta)];
M = bsxfun (@times, [-1 1; 1 1; 1 -1; -1 -1], [w h]);
tx  = x + M * v;
ty  = y + M(4:-1:1,[2 1]) * v;

if nargout <= 1
  tx = [tx ty];
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% STD_PREPARE_NEIGHBORS - prepare Fieldtrip channel neighbor structure.
%                           Only prepare the structure if necessary based
%                           on statistical options in STUDY.etc.statistics.
%                           Use the 'force' option to force preparing the
%                           matrix.
%
% Usage:
%   >> [STUDY neighbors] = std_prepare_neighbors( STUDY, ALLEEG, 'key', val)
%
% Inputs:
%   STUDY        - an EEGLAB STUDY set of loaded EEG structures
%   ALLEEG       - ALLEEG vector of one or more loaded EEG dataset structures
%
% Optional inputs:
%   'force'        - ['on'|'off'] force generating the structure irrespective
%                    of the statistics options selected. Default is 'off'.
%   'channels'     - [cell] list of channels to include in the matrix
%   'method'       - [string] method for preparing. See ft_prepare_neighbors
%   'neighbordist' - [float] max distance. See ft_prepare_neighbors
%
% Note: other ft_prepare_neighbors fields such as template, layout may
%       also be used as optional keywords.
%
% Outputs:
%   STUDY       - an EEGLAB STUDY set of loaded EEG structures
%   neighbors   - Fieldtrip channel neighbour structure
%   limostruct  - structure compatible with LIMO (need to be saved in a file)
%
% Author: Arnaud Delorme, SCCN, UCSD, 2012-
%
% See also: STATCONDFIELDTRIP

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

function [STUDY neighbors limostruct] = std_prepare_neighbors(STUDY, ALLEEG, varargin);

neighbors = [];
if nargin < 2
    return;
end

[opt addopts] = finputcheck( varargin, {  'force'    'string'  { 'on','off' }   'off';
                                          'channels' 'cell'    {}               {} }, 'std_stat', 'ignore');

if strcmpi(opt.force, 'on') || (strcmpi(STUDY.etc.statistics.fieldtrip.mcorrect, 'cluster') && ...
        strcmpi(STUDY.etc.statistics.mode, 'fieldtrip') && (strcmpi(STUDY.etc.statistics.groupstats, 'on') || strcmpi(STUDY.etc.statistics.condstats, 'on')))

    EEG = eeg_emptyset;
    EEG.chanlocs = eeg_mergelocs(ALLEEG.chanlocs);
    if isempty(EEG.chanlocs)
        disp('std_prepare_neighbors: cannot prepare channel neighbour structure because of empty channel structures');
        return;
    end

    if ~isempty(STUDY.etc.statistics.fieldtrip.channelneighbor) && isempty(addopts) && ...
        length(STUDY.etc.statistics.fieldtrip.channelneighbor) == length(EEG.chanlocs)

        disp('Using stored channel neighbour structure');
        neighbors = STUDY.etc.statistics.fieldtrip.channelneighbor;

    else

        if ~isempty(opt.channels)
            indChans = eeg_chaninds(EEG, opt.channels);
            EEG.chanlocs = EEG.chanlocs(indChans);
        end

        EEG.nbchan   = length(EEG.chanlocs);
        EEG.data     = zeros(EEG.nbchan,100,1);
        EEG.trials   = 1;
        EEG.pnts     = 100;
        EEG.xmin     = 0;
        EEG.srate    = 1;
        EEG.xmax     = 99;
        EEG = eeg_checkset(EEG);
        tmpcfg = eeglab2fieldtrip(EEG, 'preprocessing', 'none');

        % call the function that find channel neighbors
        % ---------------------------------------------
        addparams = eval( [ '{' STUDY.etc.statistics.fieldtrip.channelneighborparam '}' ]);
        for index = 1:2:length(addparams)
            tmpcfg = setfield(tmpcfg, addparams{index}, addparams{index+1});
        end
        for index = 1:2:length(addopts)
            tmpcfg = setfield(tmpcfg, addopts{index}, addopts{index+1});
        end
        warning off;
        if isfield(EEG.chanlocs, 'theta') && ~isempty(EEG.chanlocs(1).theta)
            tmpcfg = rmfield(tmpcfg, 'cfg');
            tmpcfg2 = tmpcfg;
            tmpcfg  = rmfield(tmpcfg, 'label'); % first input must not be data
            tmpcfg2 = rmfield(tmpcfg2, 'method'); % second input must not be method
            if isfield(tmpcfg, 'trialinfo')
                tmpcfg  = rmfield(tmpcfg, 'trialinfo'); % first input must not be data
            end
            % tmpcfg = rmfield(tmpcfg, 'label');
            % --> removing label seems to make ft_prepare_neighbours to crash
            cfg.neighbors = ft_prepare_neighbours(tmpcfg, tmpcfg2);
            neighbors = cfg.neighbors;
        else
            neighbors = [];
        end
        warning on;

    end

    STUDY.etc.statistics.fieldtrip.channelneighbor = neighbors;

    if nargout > 2
        limostruct.expected_chanlocs = EEG.chanlocs;

        if ~isempty(neighbors)
            allLabels = { neighbors.label };
            limostruct.channeighbstructmat = zeros(length(EEG.chanlocs));
            for iN = 1:length(neighbors)
                if ~isequal(neighbors(iN).label, limostruct.expected_chanlocs(iN).labels)
                    error('Wrong label');
                else
                    [tmp posChan] = intersect( allLabels, neighbors(iN).neighblabel);
                    limostruct.channeighbstructmat(iN,posChan) = 1;
                    limostruct.channeighbstructmat(posChan,iN) = 1;
                end
            end
        else
            limostruct.channeighbstructmat = ones(length(EEG.chanlocs));
        end
    end
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
classdef TestDrawKeylines
    %TestDrawKeylines

    properties (Constant)
        img = fullfile(mexopencv.root(),'test','building.jpg');
    end

    methods (Static)
        function test_1
            im = cv.imread(TestDrawKeylines.img, 'Color',true, 'ReduceScale',2);
            detector = cv.LSDDetector();
            klines = detector.detect(im, 'Scale',2, 'NumOctaves',1);
            out = cv.drawKeylines(im, klines);
            validateattributes(out, {class(im)}, {'size',size(im)});
        end

        function test_2
            im = cv.imread(TestDrawKeylines.img, 'Color',true, 'ReduceScale',2);
            detector = cv.LSDDetector();
            klines = detector.detect(im, 'Scale',2, 'NumOctaves',1);
            out = cv.drawKeylines(im, klines, ...
                'OutImage',im, 'Color',[255 0 0]);
            validateattributes(out, {class(im)}, {'size',size(im)});
        end

        function test_error_argnum
            try
                cv.drawKeylines();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Get the freq. and amp. parameters of a spectral peak by quandratic fitting
%
% Octave compatible
%
% Inputs
%  S        : The spectrum containing the peak to fit
%  fs       : [Hz] The sampling frequency
%  index    : The peak position where the parabola has to be fit
%  [zp]     : the zero padding factor used for the DFT computation
%             (used for the Abe & Smith corrections)
%  [wintype]: The window type used for the DFT computation
%             wintype=2 => Hann
%
% Outputs
%  freq     : [Hz] The frequency parameter
%  amp      : The amplitude parameter (linear scale)
%
% Copyright (c) 2012 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [freq amp p] = spec_fit_freq_amp(S, fs, index, zp, wintype)

    if index<=1
        freq = index;
        amp = log(abs(S(1)));
        return
    elseif index>=length(S)
        freq = index;
        amp = log(abs(S(end)));
        return
    end

    % fit a parabola in log amplitudes
    y1 = log(abs(S(index-1)));
    y2 = log(abs(S(index)));
    y3 = log(abs(S(index+1)));
    A = (y3-y2)/2 + (y1-y2)/2;
    B = -( y1 - y2 - A );
    C = y1 - A + B;
    p = [A B C];

    % max abscissa
    di = -p(2)/(2*p(1));

    % max amplitude
    logamp = p(1)*di*di + p(2)*di + p(3);

    if nargin>3 && ~isempty(zp)
        % Use Abe & Smith corrections
        if wintype==2;
            c0=0.247560; c1=0.084372; c2=-0.090608; c3=-0.055781;   % Hann
        else
            error('Unknown correction terms for the given window for Abe&Smith corrections');
        end
        delta = di;

        % Frequency correction
        ksi = c0*zp^(-2) + c1*zp^(-4);                  % (3)
        di = delta + ksi*(delta-0.5)*(delta+0.5)*delta; % (1)

        % Amplitude correction
        ita = c2*zp^(-4) + c3*zp^(-6);                  % (4)
        logamp = logamp + ita*di^2;                     % (2)
    end

    freq = index + di;
    freq = (freq-1)*(fs/length(S));

    amp = exp(logamp);
return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [ons_secs, sqpar, verbose] = tapas_physio_crop_scanphysevents_to_acq_window(ons_secs, sqpar, verbose)
% cropping of ons_secs into acquired scan/presentation session, augmenting
% sqpar by scan-timing parameters from SCANPHYSLOG-file
%
% USAGE
%   function [ons_secs, sqpar] = ...
%       tapas_physio_crop_scanphysevents_to_acq_window(ons_secs, sqpar)
%
% IN
%   ons_secs    - onsets of all physlog events in seconds
%               .spulse     = onsets of slice scan acquisition
%               .cpulse     = onsets of cardiac R-wave peaks
%               .r          = time series of respiration
%               .svolpulse  = onsets of volume scan acquisition
%               .t          = time vector of logfile rows
%
%   sqpar     - sequence timing parameters
%           .Nslices        - number of slices per volume in fMRI scan
%           .NslicesPerBeat - usually equals Nslices, unless you trigger
%                             with the heart beat
%           .TR             - repetition time in seconds
%           .Ndummies       - number of dummy volumes
%           .Nscans         - number of full volumes saved (volumes in nifti file,
%                             usually rows in your design matrix)
%            onset_slice    - slice whose scan onset determines the adjustment of the
%                             regressor timing to a particular slice for the whole volume
%
% OUT
%   ons_secs    - input ons_secs cropped to acquisition window
%                 .raw - uncropped ons_secs-structure as input into this
%                 function
%   sqpar       - augmented input, also contains
%           .maxscan        - acquired volumes during session running
%           .Nvols_paradigm - acquired volumes during paradigm running
%           .meanTR         - mean repetition time (secs) over session
%
% Lars Kasper, August 2011
% Copyright (C) 2013 Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% parameter settings
    Nscans          = sqpar.Nscans;
    Ndummies        = sqpar.Ndummies;
    NslicesPerBeat  = sqpar.NslicesPerBeat;
    Nslices         = sqpar.Nslices;

    spulse          = ons_secs.spulse;
    svolpulse       = ons_secs.svolpulse;
    cpulse          = ons_secs.cpulse;
    c               = ons_secs.c;
    r               = ons_secs.r;
    c_is_reliable = ons_secs.c_is_reliable;
    r_is_reliable = ons_secs.r_is_reliable;
    fr              = ons_secs.fr;
    t               = ons_secs.t;
    ons_secs.raw    = ons_secs;

%% cut after end of paradigm window
maxscan = Nscans + Ndummies;
tmax    = ons_secs.spulse_per_vol{maxscan}(end);

tstart  = ons_secs.spulse_per_vol{1}(1);
tend    = ons_secs.spulse_per_vol{maxscan}(end);

spulse((maxscan*Nslices+1):end) = [];

% 1st heartbeat should be right before first scan, others cut
% last heartbeat should be right after last scan; rest cut

acqwindow   = sort([find(cpulse<=tend & cpulse>=tstart); ...
                find(cpulse<tstart,1,'last'); find(cpulse>tend,1,'first')]);

if ~isempty(cpulse), cpulse      = cpulse(acqwindow); end;

% same for respiratory signal
acqwindow   = sort([find(t<=tend & t>=tstart); ...
                find(t<tstart,1,'last'); find(t>tend,1,'first')]);

if ~isempty(r),  r      = r(acqwindow); end;
if ~isempty(r_is_reliable), r_is_reliable = r_is_reliable(acqwindow); end;
if ~isempty(c_is_reliable), c_is_reliable = c_is_reliable(acqwindow); end;
if ~isempty(fr), fr     = fr(acqwindow); end;
if ~isempty(c),  c      = c(acqwindow); end;

ons_secs.t  = t(acqwindow);

% necessary vector for t1correction, all volume excitations needed
sqpar.t        = intersect(spulse,svolpulse);
if maxscan<max(size(sqpar.t))
    sqpar.t    = sqpar.t(1:maxscan*Nslices/NslicesPerBeat); %cut what's not needed
end
sqpar.maxscan = length(sqpar.t); %counts each excitation

% Counts only real repetitions of a slice
sqpar.Nvols_paradigm   = (length(sqpar.t) - ...
    Ndummies*Nslices / NslicesPerBeat)/Nslices*NslicesPerBeat;

% Mean TR of time paradigm is running, excluding dummies
lastmin1        = min([length(svolpulse)-1,sqpar.maxscan]);
sqpar.meanTR    = mean( diff(svolpulse((sqpar.Ndummies:Nslices/NslicesPerBeat:lastmin1)+1)));

formatstr = ['    maxscan (incl. dummies) = %d \n    ', ...
    'tmin (1st scan start (1st dummy))= %6.2f s\n    ',...
    'tmin (1st scan start (after dummies))= %6.2f s\n    ', ...
    'tmax = %6.2f s \n    ', ...
    'mean TR = %6.2f s\n'];

verbose = tapas_physio_log(sprintf(formatstr, sqpar.maxscan, spulse(1), ...
    spulse(1+Ndummies*Nslices), tmax, sqpar.meanTR), verbose);

%% prepare output variable

ons_secs.c          = c;
ons_secs.r          = r;
ons_secs.c_is_reliable = c_is_reliable;
ons_secs.r_is_reliable = r_is_reliable;
ons_secs.fr         = fr;
ons_secs.spulse     = spulse;
ons_secs.cpulse     = cpulse;
ons_secs.svolpulse  = svolpulse;

end

author: translationalneuromodeling

repo: tapas


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym isrow (@var{x})
%% Return true if symbolic expression is a row vector.
%%
%% Example:
%% @example
%% @group
%% h = sym([1 2 3]);
%% isrow(h)
%%   @result{} 1
%% isrow(sym(1))
%%   @result{} 1
%% isrow(h')
%%   @result{} 0
%% @end group
%% @end example
%%
%% @seealso{@@sym/iscolumn, @@sym/isvector, @@sym/isscalar}
%% @end defmethod

function r = isrow(x)

  if (nargin ~= 1)
    print_usage ();
  end

  % from Rik Wehbring's Octave function
  sz = size (x);
  r = (ndims (x) == 2 && (sz(1) == 1));

end

%!assert (isrow (sym ([1])))
%!assert (isrow (sym ([1 2 3])))
%!assert (~isrow (sym ([])))
%!assert (~isrow (sym ([1 2 3]')))
%!assert (~isrow (sym ([1 2; 3 4])))

repo: octsympy


text:
function variance = coupon_variance ( j, n )

%*****************************************************************************80
%
%% COUPON_VARIANCE returns the variance of the Coupon PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer J, the number of distinct coupons to be collected.
%    J must be between 1 and N.
%
%    Input, integer N, the number of distinct coupons.
%
%    Output, real VARIANCE, the variance of the number of
%    coupons that must be collected in order to just get J distinct kinds.
%
  if ( n < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'COUPON_VARIANCE - Fatal error!\n' );
    fprintf ( 1, '  Number of distinct coupons desired must be no more\n' );
    fprintf ( 1, '  than the total number of distinct coupons.\n' );
    error ( 'COUPON_VARIANCE - Fatal error!' );
  end

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

  return
end

repo: jburkardt-m


text:
function X = bartelsStewart(A, B, C, D, E, xsplit, ysplit)
%BARTELSSTEWART   Solution to generalized Sylvester matrix equation.
%
% Computes the solution to the Sylvester equation
%
%         AXB^T + CXD^T = E
%
% by using the Bartels--Stewart algorithm, see
%
% J. D. Gardiner, A. J. Laub, J. J. Amato, & C. B. Moler, Solution of
% the Sylvester matrix equation AXB^T+ CXD^T= E, ACM Transactions on
% Mathematical Software (TOMS), 18(2), 223-231.
%
% This Bartels--Stewart solver also takes information xsplit, ysplit so
% that if possible it decouples the even and odd modes.

% Copyright 2017 by The University of Oxford and The Chebfun2 Developers.
% See http://www.chebfun.org/ for Chebfun information.

tol = eps();

% If the RHS is zero then the solution is the zero solution (assuming
% uniqueness).
if ( norm(E) < 10*tol )
    X = zeros(size(E));
    return
end

% If the equation is even/odd in the x-direction then we can split the problem
% into two subproblems. We enforce P and S as upper triangular because they
% are up to rounding errors, and we need to do back substitution with
% them.
if ( ysplit )
    % This is equivalent to qz(full(A),full(C)), but faster.
    [P, S, Q1, Z1] = qzsplit(A, C);
    P = triu(P);
    S = triu(S);
else
    [P, S, Q1, Z1] = qz(full(A), full(C));
    P = triu(P);
    S = triu(S);
end

% If the PDE is even/odd in the y-direction then we can split (further)
% into double as many subproblems.
if ( xsplit )
    % Faster QZ when even/odd modes decouple in x-direction:
    [T, R, Q2, Z2] = qzsplit(D, B);
else
    % QZ does not take matrices in sparse format:
    [T, R, Q2, Z2] = qz(full(D), full(B));
end

% Now use the generalised Bartels--Stewart solver found in Gardiner et al.
% (1992).  The Sylvester matrix equation now contains quasi upper-triangular
% matrices and we can do a backwards substitution.

% transform the righthand side.
F = Q1*E*Q2.';

% Solution will be a m by n matrix.
m = size(A, 1);
n = size(B, 1);
Y = zeros(m, n);

% Do a backwards substitution type algorithm to construct the solution.
k=n;
PY = zeros(m);
SY = zeros(m);

% Construct columns n,n-1,...,3,2 of the transformed solution.  The first
% column is treated as special at the end.
while k > 1
    % There are two cases, either the subdiagonal contains a zero
    % T(k,k-1)=0 and then it is a backwards substitution, or T(k,k-1)~=0
    % and then we solve a 2x2 system instead.

    if ( T(k,k-1) == 0 )
        % Simple case (almost always end up here).
        rhs = F(:,k);
        if ( k < n )

            PY(:,k+1) = P*Y(:,k+1);
            SY(:,k+1) = S*Y(:,k+1);

            for jj = k+1:n
                rhs = rhs - R(k,jj)*PY(:,jj) - T(k,jj)*SY(:,jj);
            end

        end

        % find the kth column of the transformed solution.
        Y(:,k) = (R(k,k)*P + T(k,k)*S) \ rhs;

        % go to next column
        k = k-1;

    else
        % This is a straight copy from the Gardiner et al. paper, and just
        % solves for two columns at once. (works because of
        % quasi-triangular matrices.

        % Operator reduction.
        rhs1 = F(:,k-1);
        rhs2 = F(:,k);

        for jj = k+1:n
            yj = Y(:,jj);
            rhs1 = rhs1 - R(k-1,jj)*P*yj - T(k-1,jj)*S*yj;
            rhs2 = rhs2 - R(k,jj)*P*yj - T(k,jj)*S*yj;
        end

        % 2 by 2 system.
        SM = zeros(2*n);
        up = 1:n;
        down = n+1:2*n;

        SM(up,up) = R(k-1,k-1)*P + T(k-1,k-1)*S;
        SM(up,down) = R(k-1,k)*P + T(k-1,k)*S;
        SM(down,up) = R(k,k-1)*P + T(k,k-1)*S;
        SM(down,down) = R(k,k)*P + T(k,k)*S;

        % Permute the columns and rows:
        Spermuted = zeros(2*n);
        Spermuted(1:2:2*n,1:2:2*n) = SM(1:n,1:n);
        Spermuted(2:2:2*n,2:2:2*n) = SM(n+1:2*n,n+1:2*n);

        % Solve.
        UM = Spermuted \ [rhs1; rhs2];

        Y(:,k-1) = UM(up);
        Y(:,k) = UM(down);

        PY(:,k) = P*Y(:,k);
        PY(:,k-1) = P*Y(:,k-1);
        SY(:,k) = S*Y(:,k);
        SY(:,k-1) = S*Y(:,k-1);

        % We solved for two columns so go two columns further.
        k=k-2;

    end

end

if ( k == 1 )
    % Now we have just the first column to compute.
    rhs = F(:,1);
    PY(:,2) = P*Y(:,2);
    SY(:,2) = S*Y(:,2);
    for jj = 2:n
        rhs = rhs - R(1,jj)*PY(:,jj) - T(1,jj)*SY(:,jj);
    end
    Y(:,1) = (R(1,1)*P + T(1,1)*S) \ rhs;
end

% We have now computed the transformed solution so we just transform it
% back.
X = Z1*Y*Z2.';

end

function [P, S, Q1, Z1] = qzsplit(A, C)
%QZSPLIT   A faster qz factorisation for problems that decouple.
%
% This is equivalent to standard qz, except we take account of symmetry to
% reduce the computational requirements of the QZ factorisation.

% Do the QZ by splitting the problem into two subproblems.
A = full(A);

A1 = A(1:2:end,1:2:end);
C = full(C);
C1 = C(1:2:end,1:2:end);
[P1, S1, Q1, Z1] = qz(A1, C1);

A2 = A(2:2:end,2:2:end);
C2 = C(2:2:end,2:2:end);
[P2, S2, Q2, Z2] = qz(A2, C2);

[P, S, Q1, Z1] = reform(P1, P2, S1, S2, Q1, Q2, Z1, Z2);

end

function [P, S, Q, Z] = reform(P1, P2, S1, S2, Q1, Q2, Z1, Z2)
%REFORM   Recombine subproblems to form the QZ factorization.

% Initialise all the variables.
hf1 = size(P1, 1);
n = 2*hf1 - 1;
P = zeros(n);
S = zeros(n);
Q = zeros(n);
Z = zeros(n);

% Push the subproblem back together.
P(1:hf1,1:hf1) = P1;
P(hf1+1:end,hf1+1:end) = P2;

S(1:hf1,1:hf1) = S1;
S(hf1+1:end,hf1+1:end) = S2;

Q(1:hf1,1:2:end) = Q1;
Q(hf1+1:end,2:2:end) = Q2;

Z(1:2:end,1:hf1) = Z1;
Z(2:2:end,hf1+1:end) = Z2;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Multigrid Method
%
%% Subspace correction method

%% Example: MG for 1-d Poisson equation

%% MG on bisection grids: 2-D

%% MG on bisection grids: 3-D

%% FAS: MG for non-linear equation

repo: Aero-matlab


text:
function [ceq, ceqGrad] = cst_heelStrike(x0,xF,p)
% [ceq, ceqGrad] = cst_heelStrike(x0,xF,p)
%
% This function computes the heel-strike constraint as well as gradients.
%
% INPUTS:
%   t0 = time at the start of the trajectory
%   x0 = state at the start of the trajectory
%   tF = time at the end of the trajectory
%   xF = state at the end of the trajectory
%
% OUTPUTS:
%   ceq = equality constraint
%   ceqGrad = gradient of equality constraints
%

qm = xF(1:5);
qp = x0(1:5);
dqm = xF(6:10);
dqp = x0(6:10);

if nargout == 1   %numerical gradients
    [m,mi,f,fi] = autoGen_cst_heelStrike(...
        qp(1),qp(2),qp(3),qp(4),qp(5),...  %Angles "plus" - immediately after heel strike
        qm(1),qm(2),qm(3),qm(4),qm(5),...  %Angles "minus" - immediately before heel-strike
        dqm(1),dqm(2),dqm(3),dqm(4),dqm(5),...   %Rates "minus" - immediately before heel-strike
        p.m1, p.m2, p.m3, p.m4, p.m5, p.I1, p.I2, p.I3, p.I4, p.I5, p.l1, p.l2, p.l3, p.l4, p.l5, p.c1, p.c2, p.c3, p.c4, p.c5, 0);

    M = zeros(5,5);  %Mass matrix
    F = zeros(5,1);
    M(mi) = m(:,1);
    F(fi) = f(:,1);
    dqpDyn = M\F;  %Numerically invert the mass matrix

    ceqPos = qp - qm([5;4;3;2;1]);
    ceqVel = dqpDyn - dqp;
    ceq = [ceqPos;ceqVel];

else %Analytic gradients

    [m,mi,f,fi,mz,mzi,mzd,fz,fzi,fzd] = autoGen_cst_heelStrike(...
        qp(1),qp(2),qp(3),qp(4),qp(5),...  %Angles "plus" - immediately after heel strike
        qm(1),qm(2),qm(3),qm(4),qm(5),...  %Angles "minus" - immediately before heel-strike
        dqm(1),dqm(2),dqm(3),dqm(4),dqm(5),...   %Rates "minus" - immediately before heel-strike
        p.m1, p.m2, p.m3, p.m4, p.m5, p.I1, p.I2, p.I3, p.I4, p.I5, p.l1, p.l2, p.l3, p.l4, p.l5, p.c1, p.c2, p.c3, p.c4, p.c5, 0);

    M = zeros(5,5);  %Mass matrix
    F = zeros(5,1);
    M(mi) = m(:,1);
    F(fi) = f(:,1);
    dqpDyn = M\F;  %Numerically invert the mass matrix

    Mz = zeros(mzd);
    Fz = zeros(fzd);
    Mz(mzi) = mz;
    Fz(fzi) = fz;

    nz = 22;   %Number of dimensions for gradients [t0;x0;tF;xF]
    dqpDynGrad = zeros(5,nz);
    for j=1:nz
        dqpDynGrad(:,j) = M\( -Mz(:,:,j)*dqpDyn + Fz(:,:,j) );  % Derivative of a matrix inverse http://www.atmos.washington.edu/~dennis/MatrixCalculus.pdf
    end

    dqpGrad = zeros(5,nz);
    qpGrad = zeros(5,nz);
    qmGrad = zeros(5,nz);
    for i=1:5
        qpGrad(i,1+i) = 1;
        dqpGrad(i,6+i) = 1;
        qmGrad(i,12+i) = 1;
    end

    ceqPos = qp - qm([5;4;3;2;1]);
    ceqVel = dqpDyn - dqp;
    ceq = [ceqPos;ceqVel];

    ceqPosGrad = qpGrad - qmGrad([5;4;3;2;1],:);
    ceqVelGrad = dqpDynGrad - dqpGrad;
    ceqGrad = [ceqPosGrad;ceqVelGrad];

end

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function [blockedIrrRxns, sol, tol] = findBlockedIrrRxns(model, tol, varargin)
% find all blocked irreversible reactions by solving one single LP problem:
%      min sum(z_pos + z_neg)
%      s.t.   Sv = 0
%             lb <= v <= ub
%             v_j + tol*z_pos_j >= tol   for each reaction j with lb_j >= 0
%             v_j - tol*z_neg_j <= -tol  for each reaction j with ub_j <= 0
%
% USAGE:
%    [blockedIrrRxns, fluxes] = findBlockedIrrRxns(model, tol, parameters)
%
% INPUT:
%    model:           COBRA model
%
% OPTIONAL INPUTL
%    tol:             tolerance for zeros (default feasTol*10, use default if the input is smaller)
%    parameters:      COBRA and solver-specific parameters, as a input structure or parameter/value inputs
%
% OUTPUTS:
%    blockedIrrRxns:  cell array of blocked irreversible reactions
%    sol:             solution structure from solveCobraLP
%    tol:             tolerance for zeros used (might be different from the input)

if nargin < 2
    tol = 0;
end

[~, cobraParams, solverVarargin] = parseCobraVarargin(varargin, {}, {}, {}, {'LP'});
% if tol < feasTol * 10, v, z_pos, z_neg = 0 may be feasible due to tolerance
tol = max([tol, cobraParams.LP.feasTol * 10]);

LP = struct();
if 0
    rxnFwdOnly = find(model.lb >= 0);
    rxnRevOnly = find(model.ub <= 0);
else
    %avoid reactions bound at zero
    rxnFwdOnly = find(model.lb >= 0 & model.ub>0);
    rxnRevOnly = find(model.ub <= 0 & model.lb<0);
end
[nF, nR] = deal(numel(rxnFwdOnly), numel(rxnRevOnly));
[m, n] = size(model.S);
LP.A = [model.S,                        sparse(m, nF + nR); ...  % Sv = 0
    sparse(1:nF, rxnFwdOnly, 1, nF, n), sparse(1:nF, 1:nF, tol, nF, nF + nR); ...  % v + tol*z_pos >= tol
    sparse(1:nR, rxnRevOnly, 1, nR, n), sparse(1:nR, (nF + 1):(nF + nR), -tol, nR, nF + nR)];  % v - tol*z_neg <= -tol
LP.b = [model.b; tol * ones(nF, 1); -tol * ones(nR, 1)];
LP.c = [zeros(n, 1); ones(nF + nR, 1)];
LP.lb = [model.lb; zeros(nF + nR, 1)];
LP.ub = [model.ub; ones(nF + nR, 1)];
LP.csense = [repmat('E', m, 1); repmat('G', nF, 1); repmat('L', nR, 1)];
LP.osense = 1;

sol = solveCobraLP(LP, solverVarargin.LP{:});

if sol.stat ~= 1
    warning('The model is infeasible')
    blockedIrrRxns = {};
    return
end

blockedIrrRxns = model.rxns((abs(sol.full(1:n)) < tol * (1 - cobraParams.LP.optTol)) & (model.lb >= 0 | model.ub <= 0));

author: opencobra

repo: cobratoolbox