text:
function [info, str] = our_vl_simplenn_display(net, varargin)
%VL_SIMPLENN_DISPLAY  Display the structure of a SimpleNN network.
%   VL_SIMPLENN_DISPLAY(NET) prints statistics about the network NET.
%
%   INFO = VL_SIMPLENN_DISPLAY(NET) returns instead a structure INFO
%   with several statistics for each layer of the network NET.
%
%   [INFO, STR] = VL_SIMPLENN_DISPLAY(...) returns also a string STR
%   with the text that would otherwise be printed.
%
%   The function accepts the following options:
%
%   `inputSize`:: auto
%      Specifies the size of the input tensor X that will be passed to
%      the network as input. This information is used in order to
%      estiamte the memory required to process the network. When this
%      option is not used, VL_SIMPLENN_DISPLAY() tires to use values
%      in the NET structure to guess the input size:
%      NET.META.INPUTSIZE and NET.META.NORMALIZATION.IMAGESIZE
%      (assuming a batch size of one image, unless otherwise specified
%      by the `batchSize` option).
%
%   `batchSize`:: []
%      Specifies the number of data points in a batch in estimating
%      the memory consumption, overriding the last dimension of
%      `inputSize`.
%
%   `maxNumColumns`:: 18
%      Maximum number of columns in a table. Wider tables are broken
%      into multiple smaller ones.
%
%   `format`:: `'ascii'`
%      One of `'ascii'`, `'latex'`, or `'csv'`.
%
%   See also: VL_SIMPLENN().

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

opts.inputSize = [] ;
opts.batchSize = [] ;
opts.maxNumColumns = 18 ;
opts.format = 'ascii' ;
opts = vl_argparse(opts, varargin) ;

% determine input size, using first the option, then net.meta.inputSize,
% and eventually net.meta.normalization.imageSize, if any
if isempty(opts.inputSize)
  tmp = [] ;
  opts.inputSize = [NaN;NaN;NaN;1] ;
  if isfield(net, 'meta')
    if isfield(net.meta, 'inputSize')
      tmp =  net.meta.inputSize(:) ;
    elseif isfield(net.meta, 'normalization') && ...
        isfield(net.meta.normalization, 'imageSize')
      tmp = net.meta.normalization.imageSize ;
    end
    opts.inputSize(1:numel(tmp)) = tmp(:) ;
  end
end

if ~isempty(opts.batchSize)
  opts.inputSize(4) = opts.batchSize ;
end

fields={'layer', 'type', 'name', '-', ...
        'support', 'filtd', 'filtdil', 'nfilt', 'stride', 'pad', '-', ...
        'rfsize', 'rfoffset', 'rfstride', '-', ...
        'dsize', 'ddepth', 'dnum', '-', ...
        'xmem', 'wmem'};

% get the support, stride, and padding of the operators
for l = 1:numel(net.layers)
  ly = net.layers{l} ;
  switch ly.type
    case {'conv','conv_mask'}
      ks = max([size(ly.weights{1},1) ; size(ly.weights{1},2)],1) ;
      ks = (ks - 1) .* ly.dilate + 1 ;
      info.support(1:2,l) = ks ;
    case 'pool'
      info.support(1:2,l) = ly.pool(:) ;
    otherwise
      info.support(1:2,l) = [1;1] ;
  end
  if isfield(ly, 'stride')
    info.stride(1:2,l) = ly.stride(:) ;
  else
    info.stride(1:2,l) = 1 ;
  end
  if isfield(ly, 'pad')
    info.pad(1:4,l) = ly.pad(:) ;
  else
    info.pad(1:4,l) = 0 ;
  end

  % operator applied to the input image
  info.receptiveFieldSize(1:2,l) = 1 + ...
      sum(cumprod([[1;1], info.stride(1:2,1:l-1)],2) .* ...
          (info.support(1:2,1:l)-1),2) ;
  info.receptiveFieldOffset(1:2,l) = 1 + ...
      sum(cumprod([[1;1], info.stride(1:2,1:l-1)],2) .* ...
          ((info.support(1:2,1:l)-1)/2 - info.pad([1 3],1:l)),2) ;
  info.receptiveFieldStride = cumprod(info.stride,2) ;
end

% get the dimensions of the data
info.dataSize(1:4,1) = opts.inputSize(:) ;
for l = 1:numel(net.layers)
  ly = net.layers{l} ;
  if strcmp(ly.type, 'custom') && isfield(ly, 'getForwardSize')
    sz = ly.getForwardSize(ly, info.dataSize(:,l)) ;
    info.dataSize(:,l+1) = sz(:) ;
    continue ;
  end

  info.dataSize(1, l+1) = floor((info.dataSize(1,l) + ...
                                 sum(info.pad(1:2,l)) - ...
                                 info.support(1,l)) / info.stride(1,l)) + 1 ;
  info.dataSize(2, l+1) = floor((info.dataSize(2,l) + ...
                                 sum(info.pad(3:4,l)) - ...
                                 info.support(2,l)) / info.stride(2,l)) + 1 ;
  info.dataSize(3, l+1) = info.dataSize(3,l) ;
  info.dataSize(4, l+1) = info.dataSize(4,l) ;
  switch ly.type
    case 'conv'
      if isfield(ly, 'weights')
        f = ly.weights{1} ;
      else
        f = ly.filters ;
      end
      if size(f, 3) ~= 0
        info.dataSize(3, l+1) = size(f,4) ;
      end
    case {'loss', 'softmaxloss'}
      info.dataSize(3:4, l+1) = 1 ;
    case 'custom'
      info.dataSize(3,l+1) = NaN ;
  end
end

if nargout == 1, return ; end

% print table
table = {} ;
wmem = 0 ;
xmem = 0 ;
for wi=1:numel(fields)
  w = fields{wi} ;
  switch w
    case 'type', s = 'type' ;
    case 'stride', s = 'stride' ;
    case 'rfsize', s = 'rf size' ;
    case 'rfstride', s = 'rf stride' ;
    case 'rfoffset', s = 'rf offset' ;
    case 'dsize', s = 'data size' ;
    case 'ddepth', s = 'data depth' ;
    case 'dnum', s = 'data num' ;
    case 'nfilt', s = 'num filts' ;
    case 'filtd', s = 'filt dim' ;
    case 'filtdil', s = 'filt dilat' ;
    case 'wmem', s = 'param mem' ;
    case 'xmem', s = 'data mem' ;
    otherwise, s = char(w) ;
  end
  table{wi,1} = s ;

  % do input pseudo-layer
  for l=0:numel(net.layers)
    switch char(w)
      case '-', s='-' ;
      case 'layer', s=sprintf('%d', l) ;
      case 'dsize', s=pdims(info.dataSize(1:2,l+1)) ;
      case 'ddepth', s=sprintf('%d', info.dataSize(3,l+1)) ;
      case 'dnum', s=sprintf('%d', info.dataSize(4,l+1)) ;
      case 'xmem'
        a = prod(info.dataSize(:,l+1)) * 4 ;
        s = pmem(a) ;
        xmem = xmem + a ;
      otherwise
        if l == 0
          if strcmp(char(w),'type'), s = 'input';
          else s = 'n/a' ; end
        else
          ly=net.layers{l} ;
          switch char(w)
            case 'name'
              if isfield(ly, 'name')
                s=ly.name ;
              else
                s='' ;
              end
            case 'type'
              switch ly.type
                case 'normalize', s='norm';
                case 'pool'
                  if strcmpi(ly.method,'avg'), s='apool'; else s='mpool'; end
                case 'softmax', s='softmx' ;
                case 'softmaxloss', s='softmxl' ;
                otherwise s=ly.type ;
              end
            case 'nfilt'
              switch ly.type
                case 'conv'
                  if isfield(ly, 'weights'), a = size(ly.weights{1},4) ;
                  else, a = size(ly.filters,4) ; end
                  s=sprintf('%d',a) ;
                otherwise
                  s='n/a' ;
              end
            case 'filtd'
              switch ly.type
                case 'conv'
                  s=sprintf('%d',size(ly.weights{1},3)) ;
                otherwise
                  s='n/a' ;
              end
            case 'filtdil'
              switch ly.type
                case 'conv'
                  s=sprintf('%d',ly.dilate) ;
                otherwise
                  s='n/a' ;
              end

            case 'support'
              s = pdims(info.support(:,l)) ;
            case 'stride'
              s = pdims(info.stride(:,l)) ;
            case 'pad'
              s = pdims(info.pad(:,l)) ;
            case 'rfsize'
              s = pdims(info.receptiveFieldSize(:,l)) ;
            case 'rfoffset'
              s = pdims(info.receptiveFieldOffset(:,l)) ;
            case 'rfstride'
              s = pdims(info.receptiveFieldStride(:,l)) ;

            case 'wmem'
              a = 0 ;
              if isfield(ly, 'weights') ;
                for j=1:numel(ly.weights)
                  a = a + numel(ly.weights{j}) * 4 ;
                end
              end
              % Legacy code to be removed
              if isfield(ly, 'filters') ;
                a = a + numel(ly.filters) * 4 ;
              end
              if isfield(ly, 'biases') ;
                a = a + numel(ly.biases) * 4 ;
              end
              s = pmem(a) ;
              wmem = wmem + a ;
          end
        end
    end
    table{wi,l+2} = s ;
  end
end

str = {} ;
for i=2:opts.maxNumColumns:size(table,2)
  sel = i:min(i+opts.maxNumColumns-1,size(table,2)) ;
  str{end+1} = ptable(opts, table(:,[1 sel])) ;
end

table = {...
  'parameter memory', sprintf('%s (%.2g parameters)', pmem(wmem), wmem/4);
  'data memory', sprintf('%s (for batch size %d)', pmem(xmem), info.dataSize(4,1))} ;
str{end+1} = ptable(opts, table) ;

str = horzcat(str{:}) ;

if nargout == 0
  fprintf('%s', str) ;
  clear info str ;
end

% -------------------------------------------------------------------------
function str = ptable(opts, table)
% -------------------------------------------------------------------------
switch opts.format
  case 'ascii', str = pascii(table) ;
  case 'latex', str = platex(table) ;
  case 'csv',   str = pcsv(table) ;
end
str = horzcat(str,sprintf('\n')) ;

% -------------------------------------------------------------------------
function s = pmem(x)
% -------------------------------------------------------------------------
if isnan(x),       s = 'NaN' ;
elseif x < 1024^1, s = sprintf('%.0fB', x) ;
elseif x < 1024^2, s = sprintf('%.0fKB', x / 1024) ;
elseif x < 1024^3, s = sprintf('%.0fMB', x / 1024^2) ;
else               s = sprintf('%.0fGB', x / 1024^3) ;
end

% -------------------------------------------------------------------------
function s = pdims(x)
% -------------------------------------------------------------------------
if all(x==x(1))
  s = sprintf('%.4g', x(1)) ;
else
  s = sprintf('%.4gx', x(:)) ;
  s(end) = [] ;
end

% -------------------------------------------------------------------------
function str = pascii(table)
% -------------------------------------------------------------------------
str = {} ;
sizes = max(cellfun(@(x) numel(x), table),[],1) ;
for i=1:size(table,1)
  for j=1:size(table,2)
    s = table{i,j} ;
    fmt = sprintf('%%%ds|', sizes(j)) ;
    if isequal(s,'-'), s=repmat('-', 1, sizes(j)) ; end
    str{end+1} = sprintf(fmt, s) ;
  end
  str{end+1} = sprintf('\n') ;
end
str = horzcat(str{:}) ;

% -------------------------------------------------------------------------
function str = pcsv(table)
% -------------------------------------------------------------------------
str = {} ;
sizes = max(cellfun(@(x) numel(x), table),[],1) + 2 ;
for i=1:size(table,1)
  if isequal(table{i,1},'-'), continue ; end
  for j=1:size(table,2)
    s = table{i,j} ;
    str{end+1} = sprintf('%s,', ['"' s '"']) ;
  end
  str{end+1} = sprintf('\n') ;
end
str = horzcat(str{:}) ;

% -------------------------------------------------------------------------
function str = platex(table)
% -------------------------------------------------------------------------
str = {} ;
sizes = max(cellfun(@(x) numel(x), table),[],1) ;
str{end+1} = sprintf('\\begin{tabular}{%s}\n', repmat('c', 1, numel(sizes))) ;
for i=1:size(table,1)
  if isequal(table{i,1},'-'), str{end+1} = sprintf('\\hline\n') ; continue ; end
  for j=1:size(table,2)
    s = table{i,j} ;
    fmt = sprintf('%%%ds', sizes(j)) ;
    str{end+1} = sprintf(fmt, latexesc(s)) ;
    if j<size(table,2), str{end+1} = sprintf('&') ; end
  end
  str{end+1} = sprintf('\\\\\n') ;
end
str{end+1} = sprintf('\\end{tabular}\n') ;
str = horzcat(str{:}) ;

% -------------------------------------------------------------------------
function s = latexesc(s)
% -------------------------------------------------------------------------
s = strrep(s,'\','\\') ;
s = strrep(s,'_','\char`_') ;

% -------------------------------------------------------------------------
function [cpuMem,gpuMem] = xmem(s, cpuMem, gpuMem)
% -------------------------------------------------------------------------
if nargin <= 1
  cpuMem = 0 ;
  gpuMem = 0 ;
end
if isstruct(s)
  for f=fieldnames(s)'
    f = char(f) ;
    for i=1:numel(s)
      [cpuMem,gpuMem] = xmem(s(i).(f), cpuMem, gpuMem) ;
    end
  end
elseif iscell(s)
  for i=1:numel(s)
    [cpuMem,gpuMem] = xmem(s{i}, cpuMem, gpuMem) ;
  end
elseif isnumeric(s)
  if isa(s, 'single')
    mult = 4 ;
  else
    mult = 8 ;
  end
  if isa(s,'gpuArray')
    gpuMem = gpuMem + mult * numel(s) ;
  else
    cpuMem = cpuMem + mult * numel(s) ;
  end
end

repo: interpretableCNN


text:
function RFparam = RF_train(RFparam)
%% RFparam.D: input dimension
%% RFparam.M: desired output dimension
%% RFparam.gamma: bandwidth of the Gaussian kernel

%% actually, there is no training here. This function is just randomly setting the
%% code parameters.

RFparam.R = randn(RFparam.D,RFparam.M)*sqrt(RFparam.gamma);
RFparam.B = rand(1,RFparam.M) * 2 * pi;
%RFparam.T = zeros(1,RFparam.M);
RFparam.T = (rand(1,RFparam.M) * 2 - 1) * sqrt(2);

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function [rxnAssoc, metAssoc] = updateAnimalAnnotations(animalGEM, speciesSpecRxns, speciesSpecMets);
% updateAnimalAnnotations
%   Update animal GEM annotation by integrating annotiation of Human-GEM
%   and species-specific rxns/mets
%
%   Input:
%   animalGEM            a given animal GEM
%
%   speciesSpecRxns      a structure of species-specific reactions
%   speciesSpecMets      a structure of species-specific metabolites
%
%
%   Output:
%   rxnAssoc             a structure containing reaction annotations
%   metAssoc             a structure containing metabolite annotations
%
%   NOTE: The detail specification for speciesSpecRxns and speciesSpecMets
%          structures are described in addMetabolicNetwork function.
%
%   Usage: [rxnAssoc, metAssoc] = updateAnimalAnnotations(animalGEM, speciesSpecRxns, speciesSpecMets)
%

% handle input arguments
if nargin < 3 || nargin < 2
    error('Missing species-specific annotations!');
end

% checks

%% Load Human-GEM annotation files

% get model path
[ST, I]=dbstack('-completenames');
modelPath=fileparts(fileparts(ST(I).file));

% Human-GEM has to be included into Matlab path
rxnAnnotation=fullfile(modelPath,'model','reactions.tsv');
metAnnotation=fullfile(modelPath,'model','metabolites.tsv');
if isfile(rxnAnnotation) && isfile(metAnnotation)
    fprintf('Load Human-GEM annotation files.\n');
    rxnAssocHuman = importTsvFile(rxnAnnotation);
    metAssocHuman = importTsvFile(metAnnotation);
else
    error('ERROR: No annotation file(s) are found!');
end

%% integrate human and species-specific annotations

% append species-specific annotiaons
rxnAssoc = integrateAnnotation(rxnAssocHuman, speciesSpecRxns, 'rxns');
metAssoc = integrateAnnotation(metAssocHuman, speciesSpecMets, 'mets');

% remove non-existing elements and order by the fields in model
rxnAssoc = sortCleanArray(rxnAssoc, animalGEM.rxns, 'rxns');
metAssoc = sortCleanArray(metAssoc, animalGEM.mets, 'mets');

end

%% sub-functions
% integrate human and species-specific annotations
function joinedArray= integrateAnnotation(baseArray, uniqueArray, fieldName);

    joinedArray = baseArray;  % human annotation data
    nNewElements = length(uniqueArray.(fieldName));
    baseArrayFields = fieldnames(baseArray);
    uniqueArrayFields = fieldnames(uniqueArray);

    % combine annotation information
    for i = 1:length(baseArrayFields)
        n = baseArrayFields{i};
        if ismember(n,intersect(baseArrayFields, uniqueArrayFields))
            joinedArray.(n) = [baseArray.(n); uniqueArray.(n)];
        else
            if isnumeric(joinedArray.(n)(1))
                joinedArray.(n) = [baseArray.(n); zeros(nNewElements,1)];
            else
                joinedArray.(n) = [baseArray.(n); repmat({''},nNewElements,1)];
            end
        end
    end

end

% sort and clean annotation cell array according to model rxns/mets
function sortedArray= sortCleanArray(inputArray, inputCell, fieldName);

% model components (rxns or mets) must be included in the inputArray
if all(ismember(inputCell, inputArray.(fieldName)))
    idToRemove = setdiff(inputArray.(fieldName), inputCell);

    % identifiy human rxns/mets that are not in animal GEM
    if ~isempty(idToRemove)
        [~, indexToRemove] = ismember(idToRemove, inputArray.(fieldName));
        inputArray.(fieldName)(indexToRemove) = [];  %remove non-existing elements
    end

    % get the index for aligning annoation array with model rxns/mets
    [~, sortedIndex] = ismember(inputCell, inputArray.(fieldName));
    sortedArray = {};
    sortedArray.(fieldName) = inputArray.(fieldName)(sortedIndex);

    arrayFields = setdiff(fieldnames(inputArray), fieldName);
    for i = 1:length(arrayFields)
        n = arrayFields{i};
        if ~isempty(idToRemove)
	        inputArray.(n)(indexToRemove) = [];  %remove non-existing elements
        end
        sortedArray.(n) = inputArray.(n)(sortedIndex);
    end
else
    error('There are model components that are from unknown source.');
end

end

author: SysBioChalmers

repo: Human-GEM


function errorOccurred = errorHasOccurred(status)
% determines whether an error has occurred from a status struct OR cell array
% of status structs
    errorOccurred = false;
    if iscell(status)
        for iStatus = 1:numel(status)
            errorOccurred = errorOccurred || errorHasOccurred(status{iStatus});
        end
    else
        stages = getStagesFromStatus(status);
        for iStage = 1:numel(stages)
            thisStage = status.(stages{iStage});
            errorOccurred = errorOccurred || thisStage.error;
        end
    end
end


text:
function Show_Intensity(Seqs, id, para, options)

History = [Seqs(id).Time; Seqs(id).Mark];
M = round(options.Tmax./options.dt);
time_stamp = 0:options.dt:(M-1)*options.dt;
lambda = zeros(size(para.A, 1), M);

for m = 1:M
    lambda(:, m) = Intensity_HP(time_stamp(m), History, para);
end

set(gcf,'Color','w');

figure
for u = 1:size(para.A, 1)
    subplot(1,size(para.A, 1),u)
    ind = find(History(2,:) == u);
    hold on
    stem(History(1, ind), ones(1, length(ind)), 'Color', 'k', 'LineWidth', 2);
    plot(time_stamp, lambda(u,:), 'Color', [0.8, 0, 0], 'LineWidth', 2);
    hold off
    ylabel('Intensity, \lambda(t)')
    xlabel(['Event-occurrence time (' num2str(length(ind)) ' events total)'])
end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function B = boundaries(BW, conn, dir)
%BOUNDARIES Trace object boundaries.
%   B = BOUNDARIES(BW) traces the exterior boundaries of objects in
%   the binary image BW.  B is a P-by-1 cell array, where P is the
%   number of objects in the image. Each cell contains a Q-by-2
%   matrix, each row of which contains the row and column coordinates
%   of a boundary pixel.  Q is the number of boundary pixels for the
%   corresponding object.  Object boundaries are traced in the
%   clockwise direction.
%
%   B = BOUNDARIES(BW, CONN) specifies the connectivity to use when
%   tracing boundaries.  CONN may be either 8 or 4.  The default
%   value for CONN is 8.
%
%   B = BOUNDARIES(BW, CONN, DIR) specifies the direction used for
%   tracing boundaries.  DIR should be either 'cw' (trace boundaries
%   clockwise) or 'ccw' (trace boundaries counterclockwise).  If DIR
%   is omitted BOUNDARIES traces in the clockwise direction.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2003/11/21 14:22:07 $

if nargin < 3
   dir = 'cw';
end

if nargin < 2
   conn = 8;
end

L = bwlabel(BW, conn);

% The number of objects is the maximum value of L.  Initialize the
% cell array B so that each cell initially contains a 0-by-2 matrix.
numObjects = max(L(:));
if numObjects > 0
   B = {zeros(0, 2)};
   B = repmat(B, numObjects, 1);
else
   B = {};
end

% Pad label matrix with zeros.  This lets us write the
% boundary-following loop without worrying about going off the edge
% of the image.
Lp = padarray(L, [1 1], 0, 'both');

% Compute the linear indexing offsets to take us from a pixel to its
% neighbors.
M = size(Lp, 1);
if conn == 8
   % Order is N NE E SE S SW W NW.
   offsets = [-1, M - 1, M, M + 1, 1, -M + 1, -M, -M-1];
else
   % Order is N E S W.
   offsets = [-1, M, 1, -M];
end

% next_search_direction_lut is a lookup table.  Given the direction
% from pixel k to pixel k+1, what is the direction to start with when
% examining the neighborhood of pixel k+1?
if conn == 8
   next_search_direction_lut = [8 8 2 2 4 4 6 6];
else
   next_search_direction_lut = [4 1 2 3];
end

% next_direction_lut is a lookup table.  Given that we just looked at
% neighbor in a given direction, which neighbor do we look at next?
if conn == 8
   next_direction_lut = [2 3 4 5 6 7 8 1];
else
   next_direction_lut = [2 3 4 1];
end

% Values used for marking the starting and boundary pixels.
START    = -1;
BOUNDARY = -2;

% Initialize scratch space in which to record the boundary pixels as
% well as follow the boundary.
scratch = zeros(100, 1);

% Find candidate starting locations for boundaries.
[rr, cc] = find((Lp(2:end-1, :) > 0) & (Lp(1:end-2, :) == 0));
rr = rr + 1;

for k = 1:length(rr)
   r = rr(k);
   c = cc(k);
   if (Lp(r,c) > 0) & (Lp(r - 1, c) == 0) & isempty(B{Lp(r, c)})
      % We've found the start of the next boundary.  Compute its
      % linear offset, record which boundary it is, mark it, and
      % initialize the counter for the number of boundary pixels.
      idx = (c-1)*size(Lp, 1) + r;
      which = Lp(idx);

      scratch(1) = idx;
      Lp(idx) = START;
      numPixels = 1;
      currentPixel = idx;
      initial_departure_direction = [];

      done = 0;
      next_search_direction = 2;
      while ~done
         % Find the next boundary pixel.
         direction = next_search_direction;
         found_next_pixel = 0;
         for k = 1:length(offsets)
            neighbor = currentPixel + offsets(direction);
            if Lp(neighbor) ~= 0
               % Found the next boundary pixel.

               if (Lp(currentPixel) == START) & ...
                      isempty(initial_departure_direction)
                  % We are making the initial departure from
                  % the starting pixel.
                  initial_departure_direction = direction;

               elseif (Lp(currentPixel) == START) & ...
                      (initial_departure_direction == direction)
                  % We are about to retrace our path.
                  % That means we're done.
                  done = 1;
                  found_next_pixel = 1;
                  break;
               end

               % Take the next step along the boundary.
               next_search_direction = ...
                   next_search_direction_lut(direction);
               found_next_pixel = 1;
               numPixels = numPixels + 1;
               if numPixels > size(scratch, 1)
                  % Double the scratch space.
                  scratch(2*size(scratch, 1)) = 0;
               end
               scratch(numPixels) = neighbor;

               if Lp(neighbor) ~= START
                  Lp(neighbor) = BOUNDARY;
               end

               currentPixel = neighbor;
               break;
            end

            direction = next_direction_lut(direction);
         end

         if ~found_next_pixel
            % If there is no next neighbor, the object must just
            % have a single pixel.
            numPixels = 2;
            scratch(2) = scratch(1);
            done = 1;
         end
      end

      % Convert linear indices to row-column coordinates and save
      % in the output cell array.
      [row, col] = ind2sub(size(Lp), scratch(1:numPixels));
      B{which} = [row - 1, col - 1];
   end
end

if strcmp(dir, 'ccw')
   for k = 1:length(B)
      B{k} = B{k}(end:-1:1, :);
   end
end

author: HuangCongQing


text:
function [x_best, fit_best, x_pop, fit_pop stats]=complexmethod(fcn_name,bounds,gen_max,x_start,fit_start,fcn_opts,complex_opts)
%[x_best, fit_best, x_pop, fit_pop stats]=complexmethod(fcn_name,bounds,gen_max,x_start,fit_start,fcn_opts,complex_opts)
%
% Implements the Complex Method of Contrained Optimization, as proposed by Box (1965),
% improved by Guin (1968) and Krus (1992), and following the method in Andresson (2001).
%
% Takes inputs of:
% fcn_name: string including a function to maximize.  This function must give outputs of
%   [fitness x] and take inputs of (x,fcn_opts)
% bounds: 2 by n_params matrix of the minimum and maximum bounds for each parameter
% gen_max; number of times to calculate fitness before stopping (this includes the
%    initialization calculations if fit_start is not given
% x_start: a n_population by n_params matrix of initial parameters
% fit_start: a n_params vector of fitnesses of x_start params
% fcn_opts: a variable (or structure) that may be passed to the fcn_name function
% complex_opts: a variable including some options for the optimization (you probably won't need this)
%
% References
%
%Andersson J, "Multiobjective Optimization in Engineering Design - Application to fluid Power %Systems." Doctoral thesis, Division of Fluid and Mechanical Engineering Systems, Department %of Mechanical Engineering, Linkping University, 2001 %http://citeseer.ist.psu.edu/562279.html
%
%Box, M.J., "A new method of constrained optimization and a comparison with other method," %Computer Journal, Vol. 8, No. 1, pp. 42-52, 1965.
%
%Guin J. A., "Modification of the Complex method of constraint optimization," Computer Jornal, %vol. 10, pp. 416-417, 1968. 32(34)
%
%KRUS P., JANSON A., PALMBERG J.-O., “Optimization Based on Simulation for Design of Fluid %Power Systems, in Proceedings of ASME Winter Annual Meeting, Anaheim, USA, 1992.
%
%    Copyright Travis Wiens 2008
%

if nargin<2
	error('Insufficient arguements')
end
if nargin<3
	gen_max=1000;
elseif isempty(gen_max)
    gen_max=1000;
end

if nargin<4
	x_start=[];
end
if nargin<5
	fit_start=[];
end
if nargin<6
	fcn_opts=[];
end
if nargin<7
	complex_opts.alpha=1.3;%constant determines how far to overshoot centroid
	complex_opts.n_r=4;%constant for adjusting to repeated unimprovements
elseif isempty(complex_opts)
    complex_opts.alpha=1.3;%constant determines how far to overshoot centroid
	complex_opts.n_r=4;%constant for adjusting to repeated unimprovements
end

fcn_str=['[fitness, x]=' fcn_name '(x,fcn_opts);'];%string to determine fitness
%output of x allows function to change x if desired

fit_best=-inf;%initial best fitness

n_params=size(bounds,2);

gen=1;%initialize number of generations

alpha=complex_opts.alpha;%copy values
n_r=complex_opts.n_r;

if isempty(x_start)
	f_excess=1.5;%factor to increase population over minimum (must be >=1)
	n_pop=round((n_params+1)*f_excess);
	x_pop=ones(n_pop,1)*bounds(1,:) +(ones(n_pop,1)*(bounds(2,:)-bounds(1,:))).*rand(n_pop,n_params);
else
	x_pop=x_start;
	n_pop=size(x_pop,1);

end

%initialize fitness
if isempty(fit_start)
	fit_pop=zeros(1,n_pop);
	for i=1:n_pop
		x=x_pop(i,:);
		eval(fcn_str); %get fitness for each x
		x_pop(i,:)=x;
		fit_pop(i)=fitness;
		gen=gen+1;
	end
else
	fit_pop=fit_start;
	if numel(fit_pop)~=n_pop
		error('Incorrect size of fit_pop');
	end
end

stats.trace_fitness=zeros(1,gen_max);
while (gen<gen_max)

	[fit_best idx_best]=max(fit_pop);%find best fitness
	x_best=x_pop(idx_best,:);%find best params
	[fit_worst idx_worst]=min(fit_pop);%find worst
	x_worst=x_pop(idx_worst,:);%find best params

	x_centroid=(sum(x_pop)-x_worst)/(n_pop-1);%calc centroid of all but worst x

	x_new=x_centroid+alpha*(x_centroid-x_worst);%"mirror" worst point through centroid

	idx=find(x_new<bounds(1,:));%find params out of bounds
	x_new(idx)=bounds(1,idx);
	idx=find(x_new>bounds(2,:));
	x_new(idx)=bounds(2,idx);

	x=x_new;
	eval(fcn_str); %get fitness for new x
	x_new=x;
	fit_new=fitness;
	stats.trace_fitness(gen)=fitness;
	gen=gen+1;%increment generation

	if fit_new<fit_worst %check if the new point is still worst
		k_r=1;%reset number of iterations
		while (fit_new<fit_worst&gen<gen_max)
			epsilon=(n_r/(n_r+k_r-1))^((n_r+k_r-1)/n_r);
			k_r=k_r+1;%increment number of iterations
			R=rand;
			x_store=x_new;%store previous x_new
			x_new=(x_new+epsilon*x_centroid+(1-epsilon)*x_best)/2+(x_centroid-x_best)*(1-epsilon)*(2*R-1);

			idx=find(x_new<bounds(1,:));%find params out of bounds
			x_new(idx)=bounds(1,idx);
			idx=find(x_new>bounds(2,:));
			x_new(idx)=bounds(2,idx);

			x=x_new;
			eval(fcn_str); %get fitness for new x
			x_new=x;
			fit_new=fitness;
			stats.trace_fitness(gen)=fitness;
			gen=gen+1;%increment generation
		end
	end
	x_pop(idx_worst,:)=x_new;%replace worst params with new params
	fit_pop(idx_worst)=fit_new;
end

[fit_best idx_best]=max(fit_pop);%find best fitness
x_best=x_pop(idx_best,:);%find best params

author: Sable

repo: mcbench-benchmarks


text:
% crop_fmf(infmf,outfmf,startframe,endframe)
function crop_fmf(infmf,outfmf,startframe,endframe)

[header_size,version,nr,nc,bytes_per_chunk,nframes,data_format] = fmf_read_header(infmf);

if startframe < 1 || startframe > nframes || endframe < 0 || endframe > nframes || startframe > endframe,
    error('startframe or endframe is out of range.');
end

fin = fopen(infmf,'r');
%fout = fopen(outfmf,'w');
%buf = fread(fin,header_size,'uint8');
%fwrite(fout,buf);

fout = fmf_write_header(outfmf,header_size,version,nr,nc,bytes_per_chunk,endframe-startframe+1,data_format);

fseek(fin,header_size+(startframe-1)*bytes_per_chunk,'bof');
for f = startframe:endframe,
   if mod(f-startframe+1,100) == 0,
     % fprintf('Copying frame %d in range %d to %d\n',f,startframe,endframe);
   end
   buf = fread(fin,bytes_per_chunk,'uint8');
   fwrite(fout,buf);
end

fclose(fin);
fclose(fout);

author: kristinbranson

repo: JAABA


text:
function str = print(obj, inputSizes, varargin)
%PRINT Print information about the DagNN object
%   PRINT(OBJ) displays a summary of the functions and parameters in the network.
%   STR = PRINT(OBJ) returns the summary as a string instead of printing it.
%
%   PRINT(OBJ, INPUTSIZES) where INPUTSIZES is a cell array of the type
%   {'input1nam', input1size, 'input2name', input2size, ...} prints
%   information using the specified size for each of the listed inputs.
%
%   PRINT(___, 'OPT', VAL, ...) accepts the following options:
%
%   `All`:: false
%      Display all the information below.
%
%   `Layers`:: '*'
%      Specify which layers to print. This can be either a list of
%      indexes, a cell array of array names, or the string '*', meaning
%      all layers.
%
%   `Parameters`:: '*'
%      Specify which parameters to print, similar to the option above.
%
%   `Variables`:: []
%      Specify which variables to print, similar to the option above.
%
%   `Dependencies`:: false
%      Whether to display the dependency (geometric transformation)
%      of each variables from each input.
%
%   `Format`:: 'ascii'
%      Choose between `ascii`, `latex`, `csv`, 'digraph', and `dot`.
%      The first three format print tables; `digraph` uses the plot function
%      for a `digraph` (supported in MATLAB>=R2015b) and the last one
%      prints a graph  in `dot` format. In case of zero outputs, it
%      attmepts to compile and visualise the dot graph using `dot` command
%      and `start` (Windows), `display` (Linux) or `open` (Mac OSX) on your system.
%      In the latter case, all variables and layers are included in the
%      graph, regardless of the other parameters.
%
%   `FigurePath`:: 'tempname.pdf'
%      Sets the path where any generated `dot` figure will be saved. Currently,
%      this is useful only in combination with the format `dot`.
%      By default, a unique temporary filename is used (`tempname`
%      is replaced with a `tempname()` call). The extension specifies the
%      output format (passed to dot as a `-Text` parameter).
%      If not extension provided, PDF used by default.
%      Additionally, stores the .dot file used to generate the figure to
%      the same location.
%
%    `dotArgs`:: ''
%       Additional dot arguments. E.g. '-Gsize="7"' to generate a smaller
%       output (for a review of the network structure etc.).
%
%   `MaxNumColumns`:: 18
%      Maximum number of columns in each table.
%
%   See also: DAGNN, DAGNN.GETVARSIZES().

if nargin > 1 && ischar(inputSizes)
  % called directly with options, skipping second argument
  varargin = {inputSizes, varargin{:}} ;
  inputSizes = {} ;
end

opts.all = false ;
opts.format = 'ascii' ;
opts.figurePath = 'tempname.pdf' ;
opts.dotArgs = '';
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.layers = '*' ;
opts.parameters = [] ;
opts.variables = [] ;
if opts.all || nargin > 1
  opts.variables = '*' ;
end
if opts.all
  opts.parameters = '*' ;
end
opts.memory = true ;
opts.dependencies = opts.all ;
opts.maxNumColumns = 18 ;
opts = vl_argparse(opts, varargin) ;

if nargin == 1, inputSizes = {} ; end
varSizes = obj.getVarSizes(inputSizes) ;
paramSizes = cellfun(@size, {obj.params.value}, 'UniformOutput', false) ;
str = {''} ;

if strcmpi(opts.format, 'dot')
  str = printDot(obj, varSizes, paramSizes, opts) ;
  if nargout == 0
    displayDot(str, opts) ;
  end
  return ;
end

if strcmpi(opts.format,'digraph')
  str = printdigraph(obj, varSizes) ;
  return ;
end

if ~isempty(opts.layers)
  table = {'func', '-', 'type', 'inputs', 'outputs', 'params', 'pad', 'stride'} ;
  for l = select(obj, 'layers', opts.layers)
    layer = obj.layers(l) ;
    table{l+1,1} = layer.name ;
    table{l+1,2} = '-' ;
    table{l+1,3} = player(class(layer.block)) ;
    table{l+1,4} = strtrim(sprintf('%s ', layer.inputs{:})) ;
    table{l+1,5} = strtrim(sprintf('%s ', layer.outputs{:})) ;
    table{l+1,6} = strtrim(sprintf('%s ', layer.params{:})) ;
    if isprop(layer.block, 'pad')
      table{l+1,7} = pdims(layer.block.pad) ;
    else
      table{l+1,7} = 'n/a' ;
    end
    if isprop(layer.block, 'stride')
      table{l+1,8} = pdims(layer.block.stride) ;
    else
      table{l+1,8} = 'n/a' ;
    end
  end
  str{end+1} = printtable(opts, table') ;
  str{end+1} = sprintf('\n') ;
end

if ~isempty(opts.parameters)
  table = {'param', '-', 'dims', 'mem', 'fanout'} ;
  for v = select(obj, 'params', opts.parameters)
    table{v+1,1} = obj.params(v).name ;
    table{v+1,2} = '-' ;
    table{v+1,3} = pdims(paramSizes{v}) ;
    table{v+1,4} = pmem(prod(paramSizes{v}) * 4) ;
    table{v+1,5} = sprintf('%d',obj.params(v).fanout) ;
  end
  str{end+1} = printtable(opts, table') ;
  str{end+1} = sprintf('\n') ;
end

if ~isempty(opts.variables)
  table = {'var', '-', 'dims', 'mem', 'fanin', 'fanout'} ;
  for v = select(obj, 'vars', opts.variables)
    table{v+1,1} = obj.vars(v).name ;
    table{v+1,2} = '-' ;
    table{v+1,3} = pdims(varSizes{v}) ;
    table{v+1,4} = pmem(prod(varSizes{v}) * 4) ;
    table{v+1,5} = sprintf('%d',obj.vars(v).fanin) ;
    table{v+1,6} = sprintf('%d',obj.vars(v).fanout) ;
  end
  str{end+1} = printtable(opts, table') ;
  str{end+1} = sprintf('\n') ;
end

if opts.memory
  paramMem = sum(cellfun(@getMem, paramSizes)) ;
  varMem = sum(cellfun(@getMem, varSizes)) ;
  table = {'params', 'vars', 'total'} ;
  table{2,1} = pmem(paramMem) ;
  table{2,2} = pmem(varMem) ;
  table{2,3} = pmem(paramMem + varMem) ;
  str{end+1} = printtable(opts, table') ;
  str{end+1} = sprintf('\n') ;
end

if opts.dependencies
  % print variable to input dependencies
  inputs = obj.getInputs() ;
  rfs = obj.getVarReceptiveFields(inputs) ;
  for i = 1:size(rfs,1)
    table = {sprintf('rf in ''%s''', inputs{i}), '-', 'size', 'stride', 'offset'} ;
    for v = 1:size(rfs,2)
      table{v+1,1} = obj.vars(v).name ;
      table{v+1,2} = '-' ;
      table{v+1,3} = pdims(rfs(i,v).size) ;
      table{v+1,4} = pdims(rfs(i,v).stride) ;
      table{v+1,5} = pdims(rfs(i,v).offset) ;
    end
    str{end+1} = printtable(opts, table') ;
    str{end+1} = sprintf('\n') ;
  end
end

% finish
str = horzcat(str{:}) ;
if nargout == 0,
  fprintf('%s',str) ;
  clear str ;
end

end

% -------------------------------------------------------------------------
function str = printtable(opts, table)
% -------------------------------------------------------------------------
str = {''} ;
for i=2:opts.maxNumColumns:size(table,2)
  sel = i:min(i+opts.maxNumColumns-1,size(table,2)) ;
  str{end+1} = printtablechunk(opts, table(:, [1 sel])) ;
  str{end+1} = sprintf('\n') ;
end
str = horzcat(str{:}) ;
end

% -------------------------------------------------------------------------
function str = printtablechunk(opts, table)
% -------------------------------------------------------------------------
str = {''} ;
switch opts.format
  case 'ascii'
    sizes = max(cellfun(@(x) numel(x), table),[],1) ;
    for i=1:size(table,1)
      for j=1:size(table,2)
        s = table{i,j} ;
        fmt = sprintf('%%%ds|', sizes(j)) ;
        if isequal(s,'-'), s=repmat('-', 1, sizes(j)) ; end
        str{end+1} = sprintf(fmt, s) ;
      end
      str{end+1} = sprintf('\n') ;
    end

  case 'latex'
    sizes = max(cellfun(@(x) numel(x), table),[],1) ;
    str{end+1} = sprintf('\\begin{tabular}{%s}\n', repmat('c', 1, numel(sizes))) ;
    for i=1:size(table,1)
      if isequal(table{i,1},'-'), str{end+1} = sprintf('\\hline\n') ; continue ; end
      for j=1:size(table,2)
        s = table{i,j} ;
        fmt = sprintf('%%%ds', sizes(j)) ;
        str{end+1} = sprintf(fmt, latexesc(s)) ;
        if j<size(table,2), str{end+1} = sprintf('&') ; end
      end
      str{end+1} = sprintf('\\\\\n') ;
    end
    str{end+1}= sprintf('\\end{tabular}\n') ;

  case 'csv'
    sizes = max(cellfun(@(x) numel(x), table),[],1) + 2 ;
    for i=1:size(table,1)
      if isequal(table{i,1},'-'), continue ; end
      for j=1:size(table,2)
        s = table{i,j} ;
        fmt = sprintf('%%%ds,', sizes(j)) ;
        str{end+1} = sprintf(fmt, ['"' s '"']) ;
      end
      str{end+1} = sprintf('\n') ;
    end

  otherwise
    error('Uknown format %s', opts.format) ;
end
str = horzcat(str{:}) ;
end

% -------------------------------------------------------------------------
function s = latexesc(s)
% -------------------------------------------------------------------------
s = strrep(s,'\','\\') ;
s = strrep(s,'_','\char`_') ;
end

% -------------------------------------------------------------------------
function s = pmem(x)
% -------------------------------------------------------------------------
if isnan(x),       s = 'NaN' ;
elseif x < 1024^1, s = sprintf('%.0fB', x) ;
elseif x < 1024^2, s = sprintf('%.0fKB', x / 1024) ;
elseif x < 1024^3, s = sprintf('%.0fMB', x / 1024^2) ;
else               s = sprintf('%.0fGB', x / 1024^3) ;
end
end

% -------------------------------------------------------------------------
function s = pdims(x)
% -------------------------------------------------------------------------
if all(isnan(x))
  s = 'n/a' ;
  return ;
end
if all(x==x(1))
  s = sprintf('%.4g', x(1)) ;
else
  s = sprintf('%.4gx', x(:)) ;
  s(end) = [] ;
end
end

% -------------------------------------------------------------------------
function x = player(x)
% -------------------------------------------------------------------------
if numel(x) < 7, return ; end
if x(1:6) == 'dagnn.', x = x(7:end) ; end
end

% -------------------------------------------------------------------------
function m = getMem(sz)
% -------------------------------------------------------------------------
m = prod(sz) * 4 ;
if isnan(m), m = 0 ; end
end

% -------------------------------------------------------------------------
function sel = select(obj, type, pattern)
% -------------------------------------------------------------------------
if isnumeric(pattern)
  sel = pattern ;
else
  if isstr(pattern)
    if strcmp(pattern, '*')
      sel = 1:numel(obj.(type)) ;
      return ;
    else
      pattern = {pattern} ;
    end
  end
  sel = find(cellfun(@(x) any(strcmp(x, pattern)), {obj.(type).name})) ;
end
end

% -------------------------------------------------------------------------
function h = printdigraph(net, varSizes)
% -------------------------------------------------------------------------
if exist('digraph') ~= 2
  error('MATLAB graph support not present.');
end
s = []; t = []; w = [];
varsNames = {net.vars.name};
layerNames = {net.layers.name};
numVars = numel(varsNames);
spatSize = cellfun(@(vs) vs(1), varSizes);
spatSize(isnan(spatSize)) = 1;
varChannels = cellfun(@(vs) vs(3), varSizes);
varChannels(isnan(varChannels)) = 0;

for li = 1:numel(layerNames)
  l = net.layers(li); lidx = numVars + li;
  s = [s l.inputIndexes];
  t = [t lidx*ones(1, numel(l.inputIndexes))];
  w = [w spatSize(l.inputIndexes)];
  s = [s lidx*ones(1, numel(l.outputIndexes))];
  t = [t l.outputIndexes];
  w = [w spatSize(l.outputIndexes)];
end
nodeNames = [varsNames, layerNames];
g = digraph(s, t, w);
lw = 5*g.Edges.Weight/max([g.Edges.Weight; 5]);
h = plot(g, 'NodeLabel', nodeNames, 'LineWidth', lw);
highlight(h, numVars+1:numVars+numel(layerNames), 'MarkerSize', 8, 'Marker', 's');
highlight(h, 1:numVars, 'MarkerSize', 5, 'Marker', 's');
cmap = copper;
varNvalRel = varChannels./max(varChannels);
for vi = 1:numel(varChannels)
  highlight(h, vi, 'NodeColor', cmap(max(round(varNvalRel(vi)*64), 1),:));
end
axis off;
layout(h, 'force');
end

% -------------------------------------------------------------------------
function str = printDot(net, varSizes, paramSizes, otps)
% -------------------------------------------------------------------------
str = {} ;
str{end+1} = sprintf('digraph DagNN {\n\tfontsize=12\n') ;
font_style = 'fontsize=12 fontname="helvetica"';

for v = 1:numel(net.vars)
  label=sprintf('{{%s} | {%s | %s }}', net.vars(v).name, pdims(varSizes{v}), pmem(4*prod(varSizes{v}))) ;
  str{end+1} = sprintf('\tvar_%s [label="%s" shape=record style="solid,rounded,filled" color=cornsilk4 fillcolor=beige %s ]\n', ...
    net.vars(v).name, label, font_style) ;
end

for p = 1:numel(net.params)
  label=sprintf('{{%s} | {%s | %s }}', net.params(p).name, pdims(paramSizes{p}), pmem(4*prod(paramSizes{p}))) ;
  str{end+1} = sprintf('\tpar_%s [label="%s" shape=record style="solid,rounded,filled" color=lightsteelblue4 fillcolor=lightsteelblue %s ]\n', ...
    net.params(p).name, label, font_style) ;
end

for l = 1:numel(net.layers)
  label = sprintf('{ %s | %s }', net.layers(l).name, class(net.layers(l).block)) ;
  str{end+1} = sprintf('\t%s [label="%s" shape=record style="bold,filled" color="tomato4" fillcolor="tomato" %s ]\n', ...
    net.layers(l).name, label, font_style) ;
  for i = 1:numel(net.layers(l).inputs)
    str{end+1} = sprintf('\tvar_%s->%s [weight=10]\n', ...
      net.layers(l).inputs{i}, ...
      net.layers(l).name) ;
  end
  for o = 1:numel(net.layers(l).outputs)
    str{end+1} = sprintf('\t%s->var_%s [weight=10]\n', ...
      net.layers(l).name, ...
      net.layers(l).outputs{o}) ;
  end
  for p = 1:numel(net.layers(l).params)
    str{end+1} = sprintf('\tpar_%s->%s [weight=1]\n', ...
      net.layers(l).params{p}, ...
      net.layers(l).name) ;
  end
end

str{end+1} = sprintf('}\n') ;
str = cat(2,str{:}) ;
end

% -------------------------------------------------------------------------
function displayDot(str, opts)
% -------------------------------------------------------------------------
%mwdot = fullfile(matlabroot, 'bin', computer('arch'), 'mwdot') ;
dotPaths = {'/opt/local/bin/dot', 'dot'} ;
if ismember(computer, {'PCWIN64', 'PCWIN'})
  winPath = 'c:\Program Files (x86)';
  dpath = dir(fullfile(winPath, 'Graphviz*'));
  if ~isempty(dpath)
    dotPaths = [{fullfile(winPath, dpath.name, 'bin', 'dot.exe')}, dotPaths];
  end
end
dotExe = '' ;
for i = 1:numel(dotPaths)
  [~,~,ext] = fileparts(dotPaths{i});
  if exist(dotPaths{i},'file') && ~strcmp(ext, '.m')
    dotExe = dotPaths{i} ;
    break;
  end
end
if isempty(dotExe)
  warning('Could not genereate a figure because the `dot` utility could not be found.') ;
  return ;
end

[path, figName, ext] = fileparts(opts.figurePath) ;

if isempty(ext), ext = '.pdf' ; end
if strcmp(figName, 'tempname')
  figName = tempname();
end
in = fullfile(path, [ figName, '.dot' ]) ;
out = fullfile(path, [ figName, ext ]) ;

f = fopen(in, 'w') ; fwrite(f, str) ; fclose(f) ;

cmd = sprintf('"%s" -T%s %s -o "%s" "%s"', dotExe, ext(2:end), ...
  opts.dotArgs, out, in) ;
[status, result] = system(cmd) ;
if status ~= 0
  error('Unable to run %s\n%s', cmd, result) ;
end
if ~isempty(strtrim(result))
  fprintf('Dot output:\n%s\n', result) ;
end

%f = fopen(out,'r') ; file=fread(f, 'char=>char')' ; fclose(f) ;
switch computer
  case {'PCWIN64', 'PCWIN'}
    system(sprintf('start "" "%s"', out)) ;
  case 'MACI64'
    system(sprintf('open "%s"', out)) ;
  case 'GLNXA64'
    system(sprintf('display "%s"', out)) ;
  otherwise
    fprintf('The figure saved at "%s"\n', out) ;
end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [Sx,tau,fs] = stft_fw(t, x, opt)
% STFT  Compute the short-time Fourier transform
%
% Very closely based on Steven Schimmel's stft.m and istft.m from
% his SPHSC 503: Speech Signal Processing course at Univ. Washington.

if nargin<3, opt = struct(); end

if ~isfield(opt, 'window'), opt.window = round(length(x)/16); end
if length(opt.window) == 1, opt.window = hamming(opt.window); end
%if ~isfield(opt, 'overlap'), opt.overlap = length(opt.window)-1; end
opt.overlap = length(opt.window)-1;

window = opt.window / norm(opt.window, 2); % Unit norm

% Pre-pad signal
n = length(x);

% extract short-time segments from the input
Nwin = length(window);
xbuf = buffer(x, Nwin, opt.overlap, 'nodelay');

% apply the specified window
xbuf = diag(sparse(window)) * xbuf;

% take the FFT, and keep only 0 to pi
Sx = fft(xbuf,[],1);
Sx = Sx(1:floor(Nwin/2)+1, :);

% define time and frequency sampling vectors
tau = [0:size(Sx,2)-1]*(Nwin-opt.overlap) + (Nwin-1)/2;
fs = linspace(0,1,Nwin+1); fs = fs(1:floor(Nwin/2)+1);

% correct for sampling frequency
tdiff = median(diff(t));
tau = tau * tdiff;
fs = fs / tdiff;

author: HeLiangHIT


text:
function varargout = drawPolygon (px, varargin)
%DRAWPOLYGON Draw a polygon specified by a list of points.
%
%   drawPolygon(POLY);
%   Packs coordinates in a single N-by-2 array, with N the vertex number.
%
%   drawPolygon(PX, PY);
%   Specifies coordinates in separate arrays. Both array must be N-by-1,
%   with N the number of vertices.
%
%   drawPolygon(POLYS)
%   Packs coordinate of several polygons in a cell array. Each element of
%   the array is a Ni-by-2 double array.
%
%   drawPolygon(..., NAME, VALUE);
%   Specifies drawing options by using one or several parameter name-value
%   pairs, see the doc of plot function for details.
%
%   drawPolygon(AX, ...)
%   Specifies the axis to draw the polygon on.
%
%   H = drawPolygon(...);
%   Also return a handle to the list of line objects.
%
%   Example
%     % draw a red rectangle
%     poly = [10 10;40 10;40 30;10 30];
%     figure; drawPolygon(poly, 'r');
%     axis equal; axis([0 50 0 50]);
%
%     % Draw two squares
%     px = [10 20 20 10 NaN 30 40 40 30]';
%     py = [10 10 20 20 NaN 10 10 20 20]';
%     figure;
%     drawPolygon([px py], 'lineWidth', 2);
%     axis equal; axis([0 50 0 50]);
%
%   See also
%   polygons2d, drawPolyline

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-05-05
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

% Store hold state
state = ishold(gca);
hold on;

%% Check input
if nargin < 1
    error ('should specify at least one argument');
end

% check for empty polygons
if isempty (px)
    return
end

% extract handle of axis to draw on
ax = gca;
if isAxisHandle(px)
    ax          = px;
    px          = varargin{1};
    varargin(1) = [];
end

%% Manage cell arrays of polygons

% case of a set of polygons stored in a cell array
if iscell(px)
    h   = cellfun(@(x) drawPolygon(ax, x, varargin{:}), px, 'UniformOutput', 0);
    h   = horzcat(h{:});

else
    % Check size vs number of arguments
    if size(px, 2) == 2
        % Case of polygon specified as a N-by-2 array (most standard case)
        py = px(:, 2);
        px = px(:, 1);

    elseif size(px, 2) == 1
        % Case of polygon specified as two N-by-1 arrays with same length
        if nargin < 2 || nargin == 2 && ~isnumeric(varargin{1})
            error('Matgeom:invalid_input_arg', ...
                'Should specify either a N-by-2 array, or 2 N-by-1 vectors');
        end

        % Extract coordinates of polygon vertices
        py          = varargin{1};
        varargin(1) = [];

        if length(py) ~= length(px)
            error('Matgeom:invalid_input_arg', ...
                'X and Y coordinate arrays should have same lengths (%d,%d)', ...
                length(px), length(py))
        end

    else
        error('Matgeom:invalid_input_arg', 'Should specify a N-by-2 array');
    end

    % set default line format
    if isempty (varargin)
        varargin = {'b-'};
    end

    % Check case of polygons with holes
    if any (isnan (px(:)) )
        polygons = splitPolygons ([px py]);
        h        = drawPolygon (ax, polygons, varargin{:});

    else
        % ensure last point is the same as the first one
        px(end+1, :) = px(1,:);
        py(end+1, :) = py(1,:);

        % draw the polygon outline
        h = plot(ax, px, py, varargin{:});

    end % whether there where holes

end % whether input arg was a cell

if ~state
    hold off
end

% avoid returning argument if not required
if nargout > 0
    varargout = {h};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [H2,Z2,sLen2,iter]=FGD_H(V,W,H,Z,Lpos,Lneg,tol)

% Fast gradient descent (FGD) with Newton method for NPAF
% Copyright@Naiyang Guan and Dacheng Tao
% Arguments:
%     sLen: step length

% Calculate scaled negative gradient at H
n=size(H,2);
G=H.*(W'*(V./Z)+H*Lneg)./(sum(W)'*ones(1,n)+H*Lpos)-H;
L=Lpos-Lneg;
Z1=W*G;
a=sum(sum(L.*(G'*G)));
d=sum(sum(L.*(G'*H)))+sum(sum(Z1));
C=Z./Z1;

% Newton method
sLen=1;
for iter=1:20,
    [sum1,sum2]=SumC1(V,C,sLen);
    sLen1=sLen-(a*sLen-sum1+d)/(a+sum2);
    if abs(sLen1-sLen)<tol,
        break;
    else
        sLen=sLen1;
    end
end

% Step length to boundary of positive orthant
if min(G(:))>=0
    sLen3=Inf;
else
    C=H./G;
    sLen3=min(-C(C<0));
end

% Best step length
sLen2=max(min(sLen1,0.99*sLen3),1);
H2=H+sLen2*G;
Z2=Z+sLen2*Z1;

return;

author: andrewssobral

repo: lrslibrary


text:
function pobjs = APPfitGroundHough(x,y, imsize)
% pobjs = APPfitGroundHough(x,y, imsize)
% Fits the ground-vertical boundary with a set of polylines
%
% Input:
%   x, y: the x and y positions of ground-vertical boundaries
%   imsize: the image size
% Output:
%   pobjs(nplanes, [m b x1 x2 maxy disconnect])
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2005
% Permission granted to non-commercial enterprises for
% modification/redistribution under GNU GPL.
% Current Version: 1.0  09/30/2005

MAKE_PLOTS = 0; % whether to make displays

[x, inds] =  sort(x);
y = y(inds);

height = imsize(1);
width = imsize(2);
y = height-round(y*(height-1));
x = round(x*(width-1)+1);
%figure(1), hold off, plot(x, y, '.b'), axis equal, hold on

ty =y;
tx =x;
minpts = sqrt(width^2+height^2)/20;
distt = sqrt(width^2+height^2)/100;
mingapt = sqrt(width^2+height^2)/20;
count = 0;

pobjs = {};

% find initial line segments
while length(tx)>minpts

    % get the hough transform image for the points (tx, ty)
    gim = zeros(imsize);
    inds = round(ty+(tx-1)*height);
    gim(inds) = 1;
    gim = conv2(gim, fspecial('gaussian', 5, 1), 'same');
    theta = [0:179];
    [R, xp] = radon(gim, theta);
    theta = theta';

    % get the best line from the hough transform image
    [maxval, ind] = max(R(:));

    % get line parameters (ax + by = c)
    [rind,thind] = ind2sub(size(R),ind);
    t = -theta(thind)*pi/180;
    r = xp(rind);
    lines = [cos(t) sin(t) -r];
    cx = width/2-1;
    cy = height/2-1;
    lines(:,3) = lines(:,3) - lines(:,1)*cx - lines(:,2)*cy;
    a = lines(1, 1);
    b2 = lines(1, 2);
    c = lines(1, 3);

    % get line parameters in y = mx+b form
    m = -lines(:, 1)/lines(:, 2);
    b = -lines(:, 3)/lines(:, 2);

    % get the distance of each point from the line
    pdst = abs(a*tx+b2*ty+c)/sqrt(a^2+b2^2);
    inds = find(pdst < distt);

    % while there is a big gap take only the larger side of the gap
    sortx = sort(tx(inds));
    gapt = max((sortx(end)-sortx(1))/5, mingapt);
    endpoints(count+1, 1:2) = [sortx(1)  sortx(end)];
    %disp(['initlength: ' num2str(length(inds))])
    %disp(['gapt: ' num2str(gapt)])
    while (max(abs(sortx(1:end-1)-sortx(2:end))) > gapt)
        [maxval, maxind] = max(abs(sortx(1:end-1)-sortx(2:end)));
        if maxind > length(inds)/2
            inds = inds(1:maxind);
        else
            inds = inds((maxind+1):end);
        end
        sortx = sort(tx(inds));
        endpoints(count+1, 1:2) = [sortx(1)  sortx(end)];
        gapt = max((sortx(end)-sortx(1))/5, mingapt);
        %disp(['gapt: ' num2str(gapt)])
    end

    %disp(['length: ' num2str(length(inds))])
    % check that best line has enough points
    if length(inds) < minpts
        break;
    end

    count = count + 1;

    % plot the lines and points on each line
    if MAKE_PLOTS
        figure(1)
        color_chars = ['y'];
        ab = axis;
        col = color_chars(mod(count-1, length(color_chars))+1);
        if abs(m)<1
             %px = [ab(1:2)];
             px = endpoints(count, 1:2);
             plot(px, height-(m*px+b)+1, ['-' col], 'LineWidth', 5);
        else
             py = ab(3:4);
             plot((height-(py-b)+1)/m,  height-py+1,  ['-' col], 'LineWidth', 5);
        end
        %plot(tx(inds), height-ty(inds)+1, ['.' col])
        drawnow;
        pause(1)
    end

    % add assigned points to x and y
    % remove assigned points from tx and ty
    tx(inds) = [];
    ty(inds) = [];

    % store current line
    pobjs{count}(1, :) = [m b 0 0 a b2 c];

end

count = length(pobjs);

if ~isempty(pobjs)

    px = cell(count, 1);
    py = cell(count, 1);
    pmedx = zeros(count, 1);
    pmedy = zeros(count, 1);
    oldpinds = cell(count, 1);

    % reassign points to each object and re-estimate lines
    distt2 = distt;
    changed = 1;
    while changed
        changed = 0;
        pdsts = zeros(length(x), count);
        for p1 = 1:length(pobjs)
            a = pobjs{p1}(5);
            b = pobjs{p1}(6);
            c = pobjs{p1}(7);
            pdsts(:, p1) = abs(a*x+b*y+c)/sqrt(a^2+b^2);
            %pobjs{p1}(6:7) = [];
        end
        [mdists, inds] = min(pdsts, [], 2);
        inds2 = find(mdists<=distt2);
        for p1 = 1:length(pobjs)
            pinds = inds2(find(inds(inds2)==p1));
            pinds = pinds(find( (x(pinds)>=(endpoints(p1, 1)-mingapt)) & ...
                (x(pinds)<=(endpoints(p1, 2)+mingapt))));
            if length(setdiff(pinds, oldpinds{p1}))>2
                changed = 1;
                %disp('changed')
                oldpinds{p1} = pinds;
            end
            px{p1} = x(pinds);
            py{p1} = y(pinds);

            if length(pinds) > 0
                endpoints(p1, 1) = min(px{p1});
                endpoints(p1, 2) = max(px{p1});

                % re-estimate line params
                p = polyfit(px{p1},py{p1},1);
                pobjs{p1}(1:2) = p;

                pm = p(1);
                pb = p(2);

                a = 1;
                b = -1/pm;
                c = pb/pm;
                pobjs{p1}(5:7) = [a b c];

                pmedx(p1) = median(px{p1});
                pmedy(p1) = median(py{p1});

            else
                endpoints(p1, 1:2) = 0;
                px{p1} = [];
                py{p1} = [];
            end

        end
    end

    % remove smaller of overlapping segments
    reminds = [];
    for p1 = 1:length(pobjs)
        for p2 = (p1+1):length(pobjs)
            if (endpoints(p1, 1)-endpoints(p2, 1))*...
                    (endpoints(p1, 2)-endpoints(p2, 2)) < 0
                if MAKE_PLOTS, disp('removing overlapping segment'), end;
                if length(px{p1}) < length(px{p2})
                    reminds(end+1) = p1;
                else
                    reminds(end+1) = p2;
                end
            end
        end
        if isempty(px{p1})
            if MAKE_PLOTS, disp('removing empty segment'), end;
            reminds = union(reminds, p1);
        end
    end
    [pobjs, pmedx, pmedy, px, py] = removeIndices(reminds, pobjs, pmedx, pmedy, px, py);
    endpoints(reminds, :) = [];

    if MAKE_PLOTS
        figure(2)
        disp(length(pobjs))
        for p1 = 1:length(pobjs)
            color_chars = ['r' 'y' 'g' 'c' 'k'];
            ab = axis;
            col = color_chars(mod(p1-1, length(color_chars))+1);
            m = pobjs{p1}(1);
            b = pobjs{p1}(2);
            if abs(m)<1
                 tpx = [ab(1:2)];
                 plot(tpx, height-(m*tpx+b)+1, ['-' col], 'LineWidth', 3);
            else
                 tpy = ab(3:4);
                 plot((height-(tpy-b)+1)/m,  height-tpy+1,  ['-' col], 'LineWidth', 3);
            end
            plot(px{p1}, height-py{p1}+1, ['.' col])
            %plot(pmedx(p1), height-pmedy(p1)+1, '+k');
        end
        drawnow;
        pause(5)
    end

    count = length(pobjs);

    if length(pobjs) > 0

    % determine which objects should be merged
    domerge = zeros(count, count);
    intersectx = zeros(count, count);
    for p1 = 1:length(pobjs)
        for p2 = (p1+1):length(pobjs)
            % if slopes are not same
            if pobjs{p1}(1) ~= pobjs{p2}(1)
                xi = (pobjs{p2}(2)-pobjs{p1}(2))/(pobjs{p1}(1)-pobjs{p2}(1));
                yi = xi*pobjs{p1}(1) + pobjs{p1}(2);
                %plot(xi, height-yi+1, '*k')
                % if med x's appear on opposite sides of intersection
                if (sign(pmedx(p1)-xi)~=sign(pmedx(p2)-xi))
                    domerge(p1, p2) = 1;
                    domerge(p2, p1) = 1;
                    intersectx(p1, p2) = xi;
                    intersectx(p2, p1) = xi;
                end
            end
        end
    end
    intersectx = round(intersectx);
    if MAKE_PLOTS, disp(['num merges: ' num2str(sum(domerge(:))/2)]), end;

    % merge objects that should be merged
    reminds = [];
    for p1 = 1:length(pobjs)
        mergeinds = find(domerge(p1, :));
        if ~isempty(mergeinds)
            domerge(p1, mergeinds) = 0;
            domerge(mergeinds, p1) = 0;
            oldlength = 0;
            while length(mergeinds)~=oldlength
                oldlength = length(mergeinds);
                for p2 = mergeinds
                    mergeinds = union(mergeinds, find(domerge(p2, :)));
                    domerge(p2, mergeinds) = 0;
                    domerge(mergeinds, p2) = 0;
                end
            end
            mergeinds = union(mergeinds, p1);

            [tmp, sind] = sort(pmedx(mergeinds));

            % sort mergeinds by their x-medians
            mergeinds = mergeinds(sind);
            tpx = [];
            tpy = [];
            tobj = [];
            for i = 1:length(mergeinds)

                % if (i) and (i-1) should merge or or (i == 1)
                if (i == 1 || intersectx(mergeinds(i), mergeinds(i-1))~=0)

                    if MAKE_PLOTS, disp(['merge: ' num2str(i)]), end;
                    tobj(end+1, :) = pobjs{mergeinds(i)};

                    % set endpoints to intersection points
                    if size(tobj, 1) > 1

                        tobj(end-1, 4) = intersectx(mergeinds(i), mergeinds(i-1));
                        tobj(end, 3) = tobj(end-1, 4);

                    end

                    % add points from (i) to merged segment
                    tpx = [tpx ; px{mergeinds(i)}];
                    tpy = [tpy ; py{mergeinds(i)}];
                else
                    % effectively remove (i) from merginds
                    mergeinds(i) = mergeinds(i-1);
                end
            end
            pobjs{p1} = tobj;
            [px{p1}, inds] = sort(tpx);
            py{p1} = tpy(inds);
            mergeinds = setdiff(mergeinds, p1);

            % set indices to remove
            reminds = [reminds(:) ; mergeinds(:)];
        end
    end
    % remove merged segments
    [pobjs, px, py, pmedx, pmedy] = removeIndices(reminds, pobjs, px, py, pmedx, pmedy);
    endpoints(reminds, :) = [];

    % set endpoints (that were not found by merging)
    for i = 1:length(pobjs)
        for j = 1:size(pobjs{i}, 1)
            if pobjs{i}(j, 3) == 0
                pobjs{i}(j, 3) = min(px{i});
            end
            if pobjs{i}(j, 4) == 0
                pobjs{i}(j, 4) = max(px{i});
            end
        end
    end

    % find overlapping segments and resolve
    reminds = [];
    for p1 = 1:length(pobjs)
        for p2 = p1+1:length(pobjs)
            if isempty(intersect(reminds, [p1 p2]))
                overx1 = max([px{p1}(1) px{p2}(1)]);
                overx2 = min([px{p1}(end) px{p2}(end)]);
                if (overx2 > overx1) % then p1 and p2 overlap

                    if MAKE_PLOTS, disp('warning: removing segment'), end;
                    if length(px{p1}) > length(px{p2})
                        reminds(end+1) = p2;
                    else
                        reminds(end+1) = p1;
                    end

                end
            end
        end
    end

    [pobjs, px, py, pmedx, pmedy] = removeIndices(reminds, pobjs, px, py, pmedx, pmedy);

    % make sure all ground points belong to one segment
    xstart = zeros(length(pobjs), 1);
    xend = zeros(length(pobjs), 1);
    for p1 = 1:length(pobjs)
        xstart(p1) = pobjs{p1}(1, 3);
        xend(p1) = pobjs{p1}(end, 4);
    end
    [xstart, ind] = sort(xstart);
    xend = xend(ind);
    pobjs = pobjs(ind);
    if xstart(1) ~= x(1)
        pobjs{1}(1, 3) = x(1);
    end
    for p1 = 1:length(xstart)-1
        if xend(p1)<xstart(p1+1)
            if pmedy(p1) > pmedy(p1+1) % p1 lower than than p2
                pobjs{p1}(end, 4) = xstart(p1+1)-1;
            else
                pobjs{p1+1}(1, 3) = xend(p1)+1;
            end
        end
    end
    [xend, ind] = sort(xend);
    if xend(end) ~= x(end)
        pobjs{ind(end)}(end, 4) = x(end);
    end

    end % if at least one object

end

if isempty(pobjs)

    reminds = [];
    for i =2:length(x)
        if x(i) == x(i-1)
            reminds(end+1) = i;
        end
    end
    x(reminds) = [];
    y(reminds) = [];
    pobjs{1} = piecewise_linear_spline(x, y, 3);
    pobjs{1}(:, 5) = 0;
end

% reverse height
for i = 1:length(pobjs)
    ind = find(pobjs{i}(:, 5)~=0);
    pobjs{i}(ind, 5) = height - pobjs{i}(ind, 5) + 1;
    pobjs{i}(:, 1) = -pobjs{i}(:, 1);
    pobjs{i}(:, 2) = height- pobjs{i}(:, 2)+1;
end

%figure(1) %, hold off,
%hold on, plot(oldx, height-oldy+1, '.k');
if MAKE_PLOTS
    figure(3)
    color_chars = ['r' 'g' 'c' 'b' 'm'];
    count = 0;
    for p = 1:length(pobjs)
        rc = mod(p-1, 5)+1;
        for j = 1:size(pobjs{p}, 1)
            count = count + 1;
            m = pobjs{p}(j, 1);
            b = pobjs{p}(j, 2);
            ppx = pobjs{p}(j, 3:4);
            plot(ppx, ppx*m+b, ['-' color_chars(rc)], 'LineWidth', 5);
            if j > 1
                plot([ppx(1) ppx(1)], [ppx(1)*m+b ppx(1)*m+b-50], ...
                    ['-' color_chars(rc)], 'LineWidth', 5);
            end
        end
    end
    drawnow;
end

tp = pobjs;
pobjs = zeros(0, 6);
for i = 1:length(tp)
    for j = 1:size(tp{i},1)
        pobjs(end+1, 1:5) = tp{i}(j, 1:5);
    end
    pobjs(end, 6) = 1;
end

pobjs(:, 1) = pobjs(:, 1)*(width-1)/(height-1);
pobjs(:, 2) = (pobjs(:, 2)-1)/(height-1);
pobjs(:, 3:4) = (pobjs(:, 3:4)-1)/(width-1);
pobjs(:, 5) = (pobjs(:, 5)-1)/(height-1);

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

function varargout = removeIndices(indices, varargin)

for i = 1:length(varargin)
    varargout(i) = varargin(i);
    varargout{i}(indices) = [];
end

author: Cloud-CV

repo: object-proposals


text:
function [v, xml] = loadXML(filename)
%
% Reads a XML file and produces a struct array.
%
% filename = name of the XML file
% v = struct variable containing fields of type char or numeric
%
% filename can also be an URL containing XML

% Read file
if ~strcmp(filename(1:5), 'http:');
    [fid, message] = fopen(filename,'r');
    if fid == -1; error(message); end
    xml = fread(fid, 'uint8=>char');
    fclose(fid);
else
    filename = strrep(filename, '\', '/');
    xml = urlread(filename);
end

% Remove 'new line' characters from the chain
xml = char(xml(:)');
xml = strrep(xml, char(10), '');
xml = strrep(xml, char(13), '');

v = xml2struct(xml);

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function box_display_test04 ( )

%*****************************************************************************80
%
%% BOX_DISPLAY_TEST04 plots Smolyak Clenshaw Curtis index sets.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_DISPLAY_TEST04:\n' );
  fprintf ( 1, '  Plot Smolyak Clenshaw Curtis monomial indices.\n' );
  fprintf ( 1, '\n' );

  m = 20;
  n = 20;
  title_string = 'Log2(X-1)+Log2(Y-1) <= 1';
  box_display ( m, n, @cc_0, @cc_1, title_string );
  print ( '-dpng', 'cc1.png' );
  fprintf ( 1, '  Created "cc1.png".\n' );

  m = 20;
  n = 20;
  title_string = 'Log2(X-1)+Log2(Y-1) <= 2';
  box_display ( m, n, @cc_1, @cc_2, title_string );
  print ( '-dpng', 'cc2.png' );
  fprintf ( 1, '  Created "cc2.png".\n' );

  m = 20;
  n = 20;
  title_string = 'Log2(X-1)+Log2(Y-1) <= 3';
  box_display ( m, n, @cc_2, @cc_3, title_string );
  print ( '-dpng', 'cc3.png' );
  fprintf ( 1, '  Created "cc3.png".\n' );

  m = 20;
  n = 20;
  title_string = 'Log2(X-1)+Log2(Y-1) <= 4';
  box_display ( m, n, @cc_3, @cc_4, title_string );
  print ( '-dpng', 'cc4.png' );
  fprintf ( 1, '  Created "cc4.png".\n' );

  m = 20;
  n = 20;
  title_string = 'Log2(X-1)+Log2(Y-1) <= 5';
  box_display ( m, n, @cc_4, @cc_5, title_string );
  print ( '-dpng', 'cc5.png' );
  fprintf ( 1, '  Created "cc5.png".\n' );

  m = 20;
  n = 20;
  title_string = 'Log2(X-1)+Log2(Y-1) <= 6';
  box_display ( m, n, @cc_5, @cc_6, title_string );
  print ( '-dpng', 'cc6.png' );
  fprintf ( 1, '  Created "cc6.png".\n' );

  m = 20;
  n = 20;
  title_string = 'Log2(X-1)+Log2(Y-1) <= 7';
  box_display ( m, n, @cc_6, @cc_7, title_string );
  print ( '-dpng', 'cc7.png' );
  fprintf ( 1, '  Created "cc7.png".\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Y = acos(X)
    % Symbolic inverse cosine.

    % Convert inputs to SymExpression
    % X = SymExpression(X);

    % construct the operation string
    sstr = ['ArcCos[' X.s ']'];

    % create a new object with the evaluated string
    Y = SymExpression(sstr);
end

author: ayonga

repo: frost-dev

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


text:
function edges = crackPattern2(box, points, alpha, varargin)
%CRACKPATTERN2 Create a (bounded) crack pattern tessellation
%
%   E = crackPattern2(BOX, POINTS, ALPHA)
%   create a crack propagation pattern wit following parameters :
%   - pattern is bounded by area BOX which is a polygon.
%   - each crack originates from points given in POINTS
%   - directions of each crack is given by a [NxM] array ALPHA, where M is
%   the number of rays emanating from each seed/
%   - a crack stop when it reaches another already created crack.
%   - all cracks stop when they reach the border of the frame, given by box
%   (a serie of 4 points).
%   The result is a collection of edges, in the form [x1 y1 x2 y2].
%
%   E = crackPattern2(BOX, POINTS, ALPHA, SPEED)
%   Also specify speed of propagation of each crack.
%
%
%   See the result with :
%     figure;
%     drawEdge(E);
%
%   See also drawEdge
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 25/05/2004.
%

%   HISTORY :

if ~isempty(varargin)
    speed = varargin{1};
else
    speed = ones(size(points, 1), 1);
end

% Compute line equations for each initial crack.
% The 'Inf' at the end correspond to the position of the limit.
% If an intersection point is found with another line, but whose position
% is after this value, this means that another crack stopped it before it
% reach the intersection point.
NP = size(points, 1);
lines = zeros(0, 5);
for i=1:size(alpha, 2)
    lines = [lines; points speed.*cos(alpha(:,i)) speed.*sin(alpha(:,i)) Inf*ones(NP, 1)];
end
NL = size(lines, 1);

% initialize lines for borders, but assign a very high speed, to be sure
% borders will stop all cracks.
dx = (box([2 3 4 1],1)-box([1 2 3 4],1))*max(speed)*5;
dy = (box([2 3 4 1],2)-box([1 2 3 4],2))*max(speed)*5;

% add borders to the lines set
lines = [lines ; createLine(box, dx, dy) Inf*ones(4,1)];

edges = zeros(0, 4);

while true
    modif = 0;

    % try to update each line
	for i=1:NL

        % initialize first point of edge
        edges(i, 1:2) = lines(i, 1:2);

        % compute intersections with all other lines
        pi = intersectLines(lines(i,:), lines);

        % compute position of all intersection points on the current line
        pos = linePosition(pi, lines(i,:));

        % consider points to the right (positive position), and sort them
        indr = find(pos>1e-12 & pos~=Inf);
        [posr, indr2] = sort(pos(indr));

        % look for the closest intersection to the right
        for i2=1:length(indr2)

            % index of intersected line
            il = indr(indr2(i2));

            % position of point relative to intersected line
            pos2 = linePosition(pi(il, :), lines(il, :));

            % depending on the sign of position, tests if the line2 can
            % stop the current line, or if it was stopped before
            if pos2>0
                if pos2<abs(posr(i2)) && pos2<lines(il, 5)
                    if lines(i, 5) ~= posr(i2)
                        edges(i, 3:4) = pi(il,:);
                        lines(i, 5) = posr(i2);
                        modif = 1;
                    end
                    break;
                end
            end
        end   % end processing of right points of the line

	end % end processing of all lines

    % break the infinite loop if no more modification was made
    if ~modif
        break;
    end
end

% add edges of the surronding box.
edges = [edges ; box(1,:) box(2,:) ; box(2,:) box(3,:); ...
                 box(3,:) box(4,:) ; box(4,:) box(1,:)  ];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function gamma_inc_values_test ( )

%*****************************************************************************80
%
%% GAMMA_INC_VALUES_TEST demonstrates the use of GAMMA_INC_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAMMA_INC_VALUES_TEST:\n' );
  fprintf ( 1, '  GAMMA_INC_VALUES stores values of\n' );
  fprintf ( 1, '  the incomplete Gamma function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A            X            GAMMA_INC(A)(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, x, fx ] = gamma_inc_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %12f  %24.16f\n', a, x, fx );

  end

  return
end


text:
function transparency_change(o2, multval)
% Change the transparency of blobs in an fmridisplay object
%
% :Inputs:
%
%   **multval:**
%        multiply transparency values by this.
%
% values < 1 makes blobs more transparent, > 1 makes blobs more opaque

for m = 1:length(o2.activation_maps)

    h = o2.activation_maps{m}.blobhandles;

    for i = 1:length(h)

        d1 = get(h(i), 'AlphaData');

        d1 = d1 .* multval;

        d1(d1 > 1) = 1;

        set(h(i), 'AlphaData', d1);

    end

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function test_ft_plot_box

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_plot_box

figure
ft_plot_box([-1 1 2 3], 'facecolor', 'b');
axis([-4 4 -4 4])

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function fl = luflops (L, U)
%LUFLOPS compute the flop count for sparse LU factorization
%
%  Example:
%      fl = luflops (L,U)
%
%  Given a sparse LU factorization (L and U), return the flop count required
%  by a conventional LU factorization algorithm to compute it.   L and U can
%  be either sparse or full matrices.  L must be lower triangular and U must
%  be upper triangular.  Do not attempt to use this on the permuted L from
%  [L,U] = lu (A).  Instead, use [L,U,P] = lu (A) or [L,U,P,Q] = lu (A).
%
%  Note that there is a subtle undercount in this estimate.  Suppose A is
%  completely dense, but during LU factorization exact cancellation occurs,
%  causing some of the entries in L and U to become identically zero.  The
%  flop count returned by this routine is an undercount.  There is a simple
%  way to fix this (L = spones (L) + spones (tril (A))), but the fix is partial.
%  It can also occur that some entry in L is a "symbolic" fill-in (zero in
%  A, but a fill-in entry and thus must be computed), but numerically
%  zero.  The only way to get a reliable LU factorization would be to do a
%  purely symbolic factorization of A.  This cannot be done with
%  symbfact (A, 'col').
%
%  See NA Digest, Vol 00, #50, Tuesday, Dec. 5, 2000
%
% See also symbfact

%    Copyright 1998-2007, Timothy A. Davis

Lnz = full (sum (spones (L))) - 1 ;	% off diagonal nz in cols of L
Unz = full (sum (spones (U')))' - 1 ;	% off diagonal nz in rows of U
fl = 2*Lnz*Unz + sum (Lnz) ;


text:
% ----------------------------------------------------------------------------
% function hfssAssignPE(fid, Name, ObjectList, [infGND = false])
%
% Description :
% -------------
% This function creates the VB Script necessary to assign a PEC boundary to
% the given object(s).
%
% Parameters :
% ------------
% fid     - file identifier of the HFSS script file.
% Name    - name of the PEC boundary. This will appear under "Boundaries"
%           in HFSS
% ObjList - a cell array of objects to which the PEC boundary condition will
%           be applied.
% infGND  - (boolean) specify as true to make the PEC represent an infinite
%           ground plane (default is false).
%
%
% Note :
% ------
%
% Example :
% ---------
% fid = fopen('myantenna.vbs', 'wt');
% ...
% hfssAssignPE(fid, 'GNDplane', {'AntennaGND'}, true);
%
% ----------------------------------------------------------------------------

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

function hfssAssignPE(fid, Name, ObjectList, infGND)

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

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

% # of objects.
nObjects = length(ObjectList);

% create the necessary script.
fprintf(fid, '\n');
fprintf(fid, 'Set oModule = oDesign.GetModule("BoundarySetup")\n');
fprintf(fid, 'oModule.AssignPerfectE _\n');
fprintf(fid, 'Array("NAME:%s", _\n', Name);

% is infinite GND ?
if (infGND)
	fprintf(fid, '"InfGroundPlane:=", true, _\n');
else
	fprintf(fid, '"InfGroundPlane:=", false, _\n');
end

fprintf(fid, '"Objects:=", _\n');
fprintf(fid, 'Array(');
for iObj = 1:nObjects
	fprintf(fid, '"%s"', ObjectList{iObj});
	if (iObj ~= nObjects)
		fprintf(fid, ',');
	end
end
fprintf(fid, '))\n');

author: yuip

repo: hfss-api

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


text:
function D = diffmat(n, m)
%DIFFMAT   Differentiation matrices for ultraspherical spectral method.
%   D = DIFFMAT(N, M) returns the differentiation matrix that takes N Chebyshev
%   coefficients and returns N C^{(M)} coefficients that represent the derivative
%   of the Chebyshev series. Here, C^{(K)} is the ultraspherical polynomial basis
%   with parameter K.

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

if ( nargin == 1 )
    m = 1;
end

% Create the differentation matrix.
if ( m > 0 )
    D = spdiags((0 : n - 1)', 1, n, n);
    for s = 1:m-1
        D = spdiags(2*s*ones(n, 1), 1, n, n) * D;
    end
else
    D = speye(n);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [poseMat] = getGroundTruth (base_dir, imuFrames)
% Based on KITTI RAW DATA DEVELOPMENT KIT
%
% Outputs ground truth poses
%
% Input arguments:
% base_dir .... absolute path to sequence base directory (ends with _sync)
% sequence base directory

% load oxts data
oxts = loadOxtsliteData(base_dir);

% transform to poses
pose = convertOxtsToPose(oxts);

poseMat = zeros(4,4);

for i = 1:length(imuFrames)
    if i == 1
        firstPose = pose{imuFrames(i)};
    end
    poseMat(:,:,i) = inv(firstPose)*pose{imuFrames(i)};
end

repo: MSCKF


text:
function [yd3] = in32yd3(in3)
% Convert volume from cubic inches to cubic yards.
% Chad Greene 2012
yd3 = in3*0.000021433470508;


text: %XOR Mapping overload

author: marianux

repo: ecg-kit

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


text:
function [ filenames ] = palsar2_files( filename, check_nonstandard_filenames )
%PALSAR2_FILES Given one file in an ALOS PALSAR 2 data bundle of files, find the rest
%
% Written by: Wade Schwartzkopf, NGA/R
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if ~exist('check_nonstandard_filenames','var')
    % Should we look at other files in the directory by content, rather
    % than filename?
    check_nonstandard_filenames = false;
end

% A standard ALOS PALSAR 2 distribution will have several files all in the
% same directory.

FILETYPE_PREFIX = {'VOL-', 'LED-', 'IMG-[H|V]{2}-', 'TRL-'};
FILETYPE = cellfun(@(x) x(1:3), FILETYPE_PREFIX, 'UniformOutput', false);
RECTYPES.VOL = [1, 192, 192, 18, 18, 360];
RECTYPES.LED = [1, 11, 192, 18, 18, 720];
RECTYPES.IMG = [1, 50, 192, 18, 18, 720];
RECTYPES.TRL = [1, 63, 192, 18, 18, 720];

[pathstr, fnameonly, ext] = fileparts(filename);
fnameonly = [fnameonly ext];
dir_info = dir(pathstr);
for i = 1:numel(FILETYPE)
    filenames.(FILETYPE{i}) = {}; % Might be multiple
    % First check filename to see if its as expected.
    expected_filename = regexprep(fnameonly, ['^' strjoin(FILETYPE_PREFIX,'|')], FILETYPE_PREFIX(i));
    dir_info_type = dir_info(~cellfun(@isempty, regexp({dir_info.name}, expected_filename, 'once')));
    for j = 1:numel(dir_info_type)
        if check_filetype(fullfile(pathstr, dir_info_type(j).name), FILETYPE{i})
            filenames.(FILETYPE{i}){end+1} = fullfile(pathstr, dir_info_type(j).name);
        end
    end
    if check_nonstandard_filenames && isempty(filenames.(FILETYPE{i}))
        % If no expected filenames, check all files in current directory
        % for matching file type.
        % Warning: Could break if multiple datasets in the same
        % directory since we can't tell whether files are from same
        % dataset or not.  We assume files are either named
        % appropriately or grouped in directories appropriately, but
        % not necessarily both.
        for j = 1:numel(dir_info)
            if ~dir_info(j).isdir && ...
                    check_filetype(fullfile(pathstr, dir_info(j).name), FILETYPE{i})
                filenames.(FILETYPE{i}){end+1} = fullfile(pathstr, dir_info(j).name);
            end
        end
    end
    % Cleanup
    if isempty(filenames.(FILETYPE{i}))
        filenames = rmfield(filenames, FILETYPE{i});
    elseif ~strcmp(FILETYPE{i},'IMG') % No cell arrays needed for non-IMG types
        filenames.(FILETYPE{i}) = filenames.(FILETYPE{i}){1};
    end
end

    function [ bool ] = check_filetype( filename, type )
        % Read volume descriptor records
        fid=fopen(filename,'r','b');
        cur_types = [fread(fid,1,'uint32') ... % Record sequence number
            fread(fid,1,'uint8') ... % Record subtype code
            fread(fid,1,'uint8') ... % Record type
            fread(fid,1,'uint8') ... % Record subtype code
            fread(fid,1,'uint8') ... % Record subtype code
            fread(fid,1,'uint32')]; % Record length
        fclose(fid);

        bool = isequal(cur_types, RECTYPES.(type));
    end

end

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

author: ngageoint


text:
function [ u, it_num ] = multigrid_poisson_1d ( n, a, b, ua, ub, force, exact )

%*****************************************************************************80
%
%% MULTIGRID_POISSON_1D solves a 1D PDE using the multigrid method.
%
%  Discussion:
%
%    This routine solves a 1D boundary value problem of the form
%
%      - U''(X) = F(X) for A < X < B,
%
%    with boundary conditions U(A) = UA, U(B) = UB.
%
%    The multigrid method is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by William Hager.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    William Hager,
%    Applied Numerical Linear Algebra,
%    Prentice-Hall, 1988,
%    ISBN13: 978-0130412942,
%    LC: QA184.H33.
%
%  Parameters:
%
%    Input, integer N, the number of intervals.
%    N must be a power of 2.
%
%    Input, real A, B, the left and right endpoints of the region.
%
%    Input, real UA, UB, the left and right boundary values.
%
%    Input, function value = FORCE ( x ), the name of the function
%    which evaluates the right hand side.
%
%    Input, function value = EXACT ( x ), the name of the function
%    which evaluates the exact solution.
%
%    Output, integer IT_NUM, the number of iterations.
%
%    Output, real U(N+1), the computed solution.
%

%
%  Determine if we have enough storage.
%
  k = i4_log_2 ( n );

  if ( n ~= 2^k )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MULTIGRID_POISSON_1D - Fatal error!\n' );
    fprintf ( 1, '  N is not a power of 2.\n' );
    error ( 'MULTIGRID_POISSON_1D - Fatal error!' );
  end

  nl = n + n + k - 2;

  uu = zeros ( nl, 1 );
  r = zeros ( nl, 1 );
%
%  Initialization.
%
  it = 4;
  it_num = 0;
  tol = 0.0001;
  utol = 0.7;
  m = n;
%
%  Set the nodes.
%
  x = ( linspace ( a, b, n + 1 ) )';
%
%  Set the right hand side.
%
  r = zeros ( n + 1, 1 );
  r(1) = ua;
  r(2:n) = force ( x(2:n) ) / n / n;
  r(n+1) = ub;
%
%  L points to first entry of solution.
%  LL points to penultimate entry.
%
  l = 1;
  ll = n;
%
%  Gauss-Seidel iteration
%
  d1 = 0.0;
  j = 0;

  while ( 1 )

    d0 = d1;
    j = j + 1;
    [ uu(l:ll+1), d1 ] = gauss_seidel ( n + 1, r(l:ll+1), uu(l:ll+1) );
    it_num = it_num + 1;
%
%  Do at least 4 iterations at each level.
%
    if ( j < it )

      continue
%
%  Enough iterations, satisfactory decrease, on finest grid, exit.
%
    elseif ( d1 < tol && n == m )

      break
%
%  Enough iterations, satisfactory convergence, go finer.
%
    elseif ( d1 < tol )

      uu(l-1-n-n:l-1) = ctof ( n + 1, uu(l:ll+1), n + n + 1, uu(l-1-n-n:l-1) );

      n = n + n;
      ll = l - 2;
      l = l - 1 - n;
      j = 0;
%
%  Enough iterations, slow convergence, 2 < N, go coarser.
%
    elseif ( utol * d0 <= d1 && 2 < n )

      [ uu(l+n+1:l+n+1+(n/2)), r(l+n+1:l+n+1+(n/2)) ] = ...
        ftoc ( n + 1, uu(l:ll+1), r(l:ll+1), (n/2)+1 );

      n = n / 2;
      l = ll + 2;
      ll = ll + n + 1;
      j = 0;

    end

  end

  u(1:n+1) = uu(1:n+1);

  return
end


text:
function out = contains(str,pattern,~,ignoreCase)

if nargin == 4 && ignoreCase
  str = lower(str);
  pattern = lower(pattern);
end

out = any(strfind(str,pattern));

author: mtex-toolbox

repo: mtex

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


text:
function z = rdivide( x, y )

%   Disciplined convex programming information for RDIVIDE:
%      For DCP purposes, the RDIVIDE division operator X./Y is identical
%      to X.*(1./Y). The right-hand term must be constant and non-zero;
%      and if the left-hand term is nonlinear, the constat must be real.
%
%   Disciplined geometric programming information for RDIVIDE:
%      Terms in a left divide must have opposite log-curvature, so the
%      following products are permitted:
%         {log-convex} ./ {log-concave}  {log-concave} ./ {log-convex}
%         {log-affine} ./ {log-affine}
%      Note that log-affine expressions are both log-convex and
%      log-concave.
%
%   For vectors, matrices, and arrays, these rules are verified
%   indepdently for each element.

z = times( x, y, './' );

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

repo: REST


text:
function [node,elem] = cubehexmesh(cube,h)
%% CUBEHEXMESH uniform mesh of cube
%
% [node,elem] = cubehexmesh([x0,x1,y0,y1,z0,z1],h) generates a uniform
% cubic mesh of the cube [x0,x1]*[y0,y1]*[z0,z1] with mesh size h.
%
% Example
%  [node,elem] = cubehexmesh([0 1 0 2 0 3],0.5);
%
% See also: squarequadmesh, cubehexgradmesh
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

x0 = cube(1); x1 = cube(2);
y0 = cube(3); y1 = cube(4);
z0 = cube(5); z1 = cube(6);
[z,x,y] = ndgrid(z0:h:z1,x0:h:x1,y0:h:y1);
node = [x(:),y(:),z(:)];
% findnode3(node); view([83 10]);

%% Generate elements
ni = size(x,1); % number of rows
nj = size(x,2); % number of columns
nk = size(x,3); % number of pages
nij = ni*nj;    % number of elements in one page
N = size(node,1);
nodeidx = reshape(1:N,ni,nj,nk);
t2nidxMap = nodeidx(1:ni-1,1:nj-1,1:nk-1);
s = t2nidxMap(:);
elem = [s s+ni s+ni+nij s+nij s+1 s+ni+1 s+ni+nij+1 s+nij+1];
% s+1 s+ni*nj
%  | /
%  s - s + ni

repo: Aero-matlab


text:
function out_fem_msh(FemMat, OutputFile)
% OUT_FEM_MSH Write tetrahedral mesh to DGF file (used by GMSH)
%
% INPUT:
%    - TessMat    : Brainstorm FEM head model mesh (tetrahedral)
%    - OutputFile : Full path to output file
%
% DOCUMENTATION:
%    - http://www.ensta-paristech.fr/~kielbasi/docs/gmsh.pdf

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

% Check the format of the mesh
if (size(FemMat.Elements,2) ~= 4)
    error('This function handles only tetrahedral meshes.');
end
% Ensure end of lines are always represented by \r\n
if ispc
    EndLine = '\n';
else
    EndLine = '\r\n';
end
% Open file
[fid, message] = fopen(OutputFile, 'wt');
if (fid < 0)
    error(['Could not create file: ' message]);
end

% Write header
fprintf(fid, ['%s' EndLine], '$MeshFormat');
fprintf(fid, ['%s' EndLine], '2.2 0 8');
fprintf(fid, ['%s' EndLine], '$EndMeshFormat');
% Write nodes
nNodes = size(FemMat.Vertices,1);
fprintf(fid, ['%s' EndLine], '$Nodes');
fprintf(fid, ['%i' EndLine], nNodes);
fprintf(fid, ['%i %i %i %i' EndLine], [1:nNodes; FemMat.Vertices']);
fprintf(fid, ['%s' EndLine], '$EndNodes');
% Write elements
nElem = size(FemMat.Elements,1);
fprintf(fid, EndLine);
fprintf(fid, ['%s' EndLine], '$Elements');
fprintf(fid, ['%i' EndLine], nElem);
fprintf(fid, ['%i %i %i %i %i %i %i %i %i' EndLine],...
    [1:nElem; ...             % elm?number: Element index
     4 * ones(1,nElem); ...   % elm?type: 4 nodes per element (tetrahedral mesh)
     2 * ones(1,nElem); ...   % number?of?tags: Two tags per element (tissue index and unused value)
     FemMat.Tissue' - 1; ...  % Tag #1: Tissue index
     zeros(1,nElem); ...      % Tag #2: Unused
     FemMat.Elements']);      % Tetrahedral elements
fprintf(fid, ['%s' EndLine], '$EndElements');

% Close file
fclose(fid);

% % LATER: If we need to visualize data within gmsh
% % bloc physical Names
% fprintf(fid,'%s' EndLine],'$PhysicalNames');
% fprintf(fid,'%i ' EndLine],length(unique(newelem(:,6))));
% fprintf(fid,'%i %i %s' EndLine],'1 1 toto1');
% fprintf(fid,'%i %i %s' EndLine],'2 2 toto2');
% fprintf(fid,'%i %i %s' EndLine],'3 3 toto3');
% fprintf(fid,'%i %i %s' EndLine],'4 4 toto4');
% fprintf(fid,'%s' EndLine],'$EndPhysicalNames');
%
% % $NodeData
% % http://geuz.org/gmsh/doc/texinfo/gmsh.html#SEC62
% fprintf(fid,'%s' EndLine],'$NodeData');
% fprintf(fid,'%i ' EndLine],1); % one string tag:
% fprintf(fid,'%s' EndLine],'"A scalar view"');
% fprintf(fid,'%i' EndLine],1); %one real tag:
% fprintf(fid,'%i' EndLine],0.0); %the time value (0.0)
% fprintf(fid,'%i' EndLine],3); % three integer tags:
% fprintf(fid,'%i' EndLine],0); %the time step (0; time steps always start at 0)
% fprintf(fid,'%i' EndLine],1); % 1-component (scalar) field
% fprintf(fid,'%i' EndLine],length(newnode)); % nb associated nodal values
% fprintf(fid,'%i %i ' EndLine],[nn,Vn']');
% fprintf(fid,'%s' EndLine],'$End$NodeData');

author: brainstorm-tools


text:
function CavitySystemID(DataFileName)

% building the Koopman linear system for cavity flow
% from the data in the file named "DataFileName"
disp('loading EDMD data ...')

% reminder
if ~strcmp(DataFileName,'Cavity_data_4EDMD_0.mat')
        warning(['if yu are using your own data', ...
        'make sure to remove the snapshots that contain NaN or Inf']);
end

load(DataFileName)

% the results are better if we center the data around point close to action
% in the state space
SUBTRACT_MEAN = 1;
if(SUBTRACT_MEAN)
    load('CavityStateLibrary.mat','FixedPoint_Re10k')
    x_mean = FixedPoint_Re10k;
    X = bsxfun(@minus,X,x_mean);
    Y = bsxfun(@minus,Y,x_mean);
end

% the following sensor locations generate the plots in the paper
load('CavitySensorLocations.mat','pindex')

% construct models for various number (k) of spatial measurements
for k = [2,5,10,40,50,100]
    tic

    Grid= CavityGridOperators( SimPar.N );
    Grid.Interior_index = setdiff(1:length(Grid.xx),Grid.Walls);

    sampind = pindex(1:k);
    % sampind is the index of the measurement point on the Chebyshev grid
    % feel free to replace sampind with some random numbers between 1 and
    % 47^2,

    CollectOutput = @(x) x(sampind,:);

    % location in the cavity domain
    figure(311)
    plot(CollectOutput(Grid.xx),CollectOutput(Grid.yy),'x')
    axis square
    axis([-1,1,-1,1])
    title(['location of sensors in cavity, k=',num2str(k)])
    drawnow
    XX = CollectOutput(X);
    YY = CollectOutput(Y);

    % embed the data from each trajectory
    nd = 5; % delay order
    tic
    Xtilde =[]; Ytilde=[]; Utilde=[];

    SimLength=unique(SimLengthVar);
    if length(SimLength)>1
        error('trajectories of different lengtnhs .. adopt')
    end
    for i=1:length(SimLengthVar)

        % consider adopting to different trajectory lengths

        TrajIndex = (i-1)*SimLength + (1:SimLength);    % indexes of each trajectory

        % embed the data
        % this gives upside down - think if ordering matters
        Xe = DelayEmbed(XX(:,TrajIndex) , nd    );
        Ye = DelayEmbed(YY(:,TrajIndex) , nd    );
        Ue = DelayEmbed( U(:,TrajIndex) , nd - 1);

        % add the delay-embedded input as an observable
        Xtilde = [Xtilde,[Xe;Ue(:,1:end-1)]];
        Ytilde = [Ytilde,[Ye;Ue(:,2:end  )]];

        UIndex = (i-1)*SimLength + (nd:SimLength);    % indexes of each trajectory
        Utilde = [Utilde, U(:,UIndex) ];

    end

    % add l-2 norm of recent output and constant as extra observables
    KE_embed = @(y) sum(abs(y((nd-1)*k+1:nd*k,:)).^2,1);

    % now lift the data
    X_lift = [Xtilde;KE_embed(Xtilde);ones(1,size(Xtilde,2))];
    Y_lift = [Ytilde;KE_embed(Ytilde);ones(1,size(Ytilde,2))];
    U_lift = Utilde;

    % compute the Koopman linear system
    % Via EDMD in the normal form
    W = [X_lift;U_lift]*[X_lift;U_lift]';
    V = Y_lift*[X_lift;U_lift]';
    M = V*pinv(W);

    Nlift = size(X_lift,1);
    A = M(:,1:Nlift);
    B = M(:,Nlift+1:end);
    C = zeros(k,size(A,1));
    C(:,(nd-1)*k+1:nd*k)=eye(k);

    % define a lift function to be used in control:
    % first a simple dely-embed that lines up the last nd columns of given
    % data matrix
    Delay_Embed = @(X,nd) reshape(X(:,end-nd+1:end),nd*size(X,1),1);

    % then a function that forms zeta given the state and input matrices
    CreateZeta = @(X,U) [Delay_Embed(CollectOutput(X(:,end-nd+1:end)),nd)...
                            ;Delay_Embed(U(:,1:end-1),nd-1)];

    % then the second lifting
    % the last function gets the actual state and input and returns the
    % Koopman linear state
    LiftFun = @(X,U) [CreateZeta(X,U); sum(CollectOutput(X(:,end)).^2);1];

    modelfilename = ['KoopmanLinSys_Re13_k',num2str(k),'_random'];
    save(modelfilename,'A','B','C','nd','CollectOutput','KE_embed','x_mean','LiftFun','CreateZeta','Delay_Embed')

    fprintf('Koopman model of cavity created using %d measurements \n',k)
    toc
end

author: arbabiha


text:
function c = correlation_rational_quadratic ( n, rho, rho0 )

%*****************************************************************************80
%
%% CORRELATION_RATIONAL_QUADRATIC evaluates the rational quadratic correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Petter Abrahamsen,
%    A Review of Gaussian Random Fields and Correlation Functions,
%    Norwegian Computing Center, 1997.
%
%  Parameters:
%
%    Input, integer N, the number of arguments.
%
%    Input, real RHO(N,1), the arguments.
%
%    Input, real RHO0, the correlation length.
%
%    Output, real C(N,1), the correlations.
%
  rho = rho ( : );

  rhohat = rho / rho0;

  c = 1.0 ./ ( 1.0 + rhohat.^2 );

  return
end


text:
% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: 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.
%
% RTB 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.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function q = interp(Q1, Q2, r)
%Quaternion.interp Interpolate rotations expressed by quaternion objects
%
% QI = Q1.interp(Q2, R) is a unit-quaternion that interpolates between Q1 for R=0
% to Q2 for R=1. This is a spherical linear interpolation (slerp) that can be
% interpretted as interpolation along a great circle arc on a sphere.
%
% If R is a vector QI is a vector of quaternions, each element
% corresponding to sequential elements of R.
%
% Notes:
% - the value of r is clipped to the interval 0 to 1
%
% See also ctraj, Quaternion.scale.

    q1 = double(Q1);
    q2 = double(Q2);

    theta = acos(q1*q2');
    count = 1;

    % clip values of r
    r(r<0) = 0;
    r(r>1) = 1;

    if length(r) == 1
        if theta == 0
            q = Q1;
        else
            q = Quaternion( (sin((1-r)*theta) * q1 + sin(r*theta) * q2) / sin(theta) );
        end
    else
        for R=r(:)'
            if theta == 0
                qq = Q1;
            else
                qq = Quaternion( (sin((1-R)*theta) * q1 + sin(R*theta) * q2) / sin(theta) );
            end
            q(count) = qq;
            count = count + 1;
        end
    end
end

author: petercorke

repo: robotics-toolbox-matlab