text:
function [net,stats] = cnn_train_autonn(net, imdb, getBatch, varargin)
%CNN_TRAIN_AUTONN Demonstrates training a CNN using the AutoNN wrapper
%   CNN_TRAIN_AUTONN is similar to CNN_TRAIN, but works with the AutoNN
%   wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-18 Andrea Vedaldi and Joao F. Henriques.
% 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).
addpath(fullfile(vl_rootnn, 'examples'));

warning('AutoNN:cnn_train_autonn', ['The function cnn_train_autonn is ' ...
  'deprecated. It will still be supported, however switching to the new, ' ...
  'more flexible training classes (Solver, Dataset, Stats) is recommended. ' ...
  'See the examples for more information. This warning will now be disabled.']) ;
warning('off', 'AutoNN:cnn_train_autonn') ;

% this is needed to harmonize the behavior of two versions of vl_nnloss:
% the legacy behavior which *sums* the loss over the batch, and the new
% behavior that takes the *average* over the batch.
% first, detect if the new behavior ('normalise' option) is present.
old = false ;
try
  vl_nnloss([], [], 'normalise', true)  ;
catch  % unrecognized option, must be the old vl_nnloss
  old = true ;
end
if old
  warning('AutoNN:normalizedLoss', ['The most recent version of ' ...
   'vl_nnloss normalizes the loss by the batch size. The current version ' ...
   'does not. A workaround is being used, but consider updating MatConvNet.']) ;
end
% consider the loss normalized (workaround in Layer.vl_nnloss). this can
% be disabled if detrimental (e.g. custom loss without normalization).
opts.normalizedLoss = true ;

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;

opts.solver = [] ;  % Empty array means use the default SGD solver
[opts, varargin] = vl_argparse(opts, varargin) ;
if ~isempty(opts.solver)
  assert(isa(opts.solver, 'function_handle') && nargout(opts.solver) == 2,...
    'Invalid solver; expected a function handle with two outputs.') ;
  % Call without input arguments, to get default options
  opts.solverOpts = opts.solver() ;
end

opts.onNanInf = [] ;  % Behavior when variables with NaN/Inf are found: 'error', 'warning' or 'debug' (default: nothing).

opts.momentum = 0.9 ;
opts.saveSolverState = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

if isa(net, 'dagnn.DagNN')
  opts.derOutputs = {'objective', 1} ;
else
  opts.derOutputs = 1 ;
end
opts.stats = 'auto' ;  % list of layers to aggregate stats from (e.g. loss, error), or 'auto' for automatic
opts.extractStatsFn = [] ;
opts.plotStatistics = true;
opts.plotDiagnostics = false ;
opts.postEpochFn = [] ;  % postEpochFn(net,params,state) called after each epoch; can return a new learning rate, 0 to stop, [] for no change
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isscalar(opts.train) && isnumeric(opts.train) && isnan(opts.train)
  opts.train = [] ;
end
if isscalar(opts.val) && isnumeric(opts.val) && isnan(opts.val)
  opts.val = [] ;
end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  if isempty(opts.derOutputs)
    error('DEROUTPUTS must be specified when training.\n') ;
  end
end

% find loss layers now instead of at every extractStatsFn call
if isa(net, 'dagnn.DagNN')
  if isempty(opts.extractStatsFn)
    opts.extractStatsFn = @extractStats ;
  end

  if isequal(opts.stats, 'auto')  % all losses
    sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
  else
    sel = zeros(size(opts.stats)) ;
    for i = 1:numel(sel)  % find layers with the given names
      sel(i) = net.getLayerIndex(opts.stats{i}) ;
    end
  end
else  % AutoNN
  if isempty(opts.extractStatsFn)
    opts.extractStatsFn = @extractStatsAutoNN ;
  end

  if isequal(opts.stats, 'auto')  % all losses
    sel = find(cellfun(@(f) isequal(f, @vl_nnloss) || ...
      isequal(f, @vl_nnsoftmaxloss), {net.forward.func})) ;
  else
    sel = zeros(size(opts.stats)) ;
    for i = 1:numel(sel)  % find layers with the given names
      sel(i) = find(strcmp({net.forward.name}, opts.stats{i})) ;
    end
  end
end
fn = opts.extractStatsFn ;
normalizedLoss = opts.normalizedLoss ;
opts.extractStatsFn = @(stats, net, batchSize) fn(stats, net, sel, batchSize, normalizedLoss) ;

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
  prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.val = opts.val(randperm(numel(opts.val))) ;
  params.imdb = imdb ;
  params.getBatch = getBatch ;

  if numel(opts.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train') ;
    [net, state] = processEpoch(net, state, params, 'val') ;
    if ~evaluateMode
      saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    spmd
      [net, state] = processEpoch(net, state, params, 'train') ;
      [net, state] = processEpoch(net, state, params, 'val') ;
      if labindex == 1 && ~evaluateMode
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  if ~isempty(opts.postEpochFn)
    state.stats = lastStats ;  % allow the func. to read accumulated stats
    if nargout(opts.postEpochFn) <= 0
      opts.postEpochFn(net, params, state) ;
    else  % returned an updated learning rate
      if nargout(opts.postEpochFn) == 1
        lr = opts.postEpochFn(net, params, state) ;
      else  % returned updated/custom stats
        [lr, lastStats] = opts.postEpochFn(net, params, state) ;
      end
      if ~isempty(lr), opts.learningRate = lr; end
      if opts.learningRate == 0, break; end
    end
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  saveStats(modelPath(epoch), stats) ;

  if opts.plotStatistics
    switchFigure(1) ; clf ;
    plots = setdiff(...
      cat(2,...
      fieldnames(stats.train)', ...
      fieldnames(stats.val)'), {'num', 'time'}) ;
    for p = plots
      p = char(p) ;
      values = zeros(0, epoch) ;
      leg = {} ;
      for f = {'train', 'val'}
        f = char(f) ;
        if isfield(stats.(f), p)
          tmp = [stats.(f).(p)] ;
          values(end+1,:) = tmp(1,:)' ;
          leg{end+1} = f ;
        end
      end
      subplot(1,numel(plots),find(strcmp(p,plots))) ;
      plot(1:epoch, values','o-') ;
      xlabel('epoch') ;
      title(p) ;
      legend(leg{:}) ;
      grid on ;
    end
    drawnow ;
    print(1, modelFigPath, '-dpdf') ;
  end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.solverState)
  state.solverState = cell(1, numel(net.params)) ;
  state.solverState(:) = {0} ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net.move('gpu') ;
  for i = 1:numel(state.solverState)
    if ~isa(net, 'dagnn.DagNN') && ~net.isGpuVar(net.params(i).var)
      continue
    end
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gpuArray(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gpuArray, s, 'UniformOutput', false) ;
    end
  end
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;

  if ~isa(net, 'dagnn.DagNN')
    % some layers have implicit Param sizes: the Params are scalar, until
    % the first derivative computation, when they're expanded to full size.
    % this allows initializing layers without knowing their input/output
    % sizes, making them less verbose. e.g. batch-norm: y = vl_nnbnorm(x);
    % for multi-GPU we need their sizes, so compute their derivatives once.
    subset = params.(mode) ;
    inputs = params.getBatch(params.imdb, subset(1)) ;  % one sample

    net.eval(inputs, 'normal', params.derOutputs) ;

    paramDer = net.getDer([net.params.var]) ;
    if ~iscell(paramDer), paramDer = {paramDer} ; end

    net.setParameterServer(parserv, paramDer) ;
  else
    net.setParameterServer(parserv) ;
  end
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

num = 0 ;
epoch = params.epoch ;
subset = params.(mode) ;
adjustTime = 0 ;

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

start = tic ;
for t=1:params.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d:', mode, epoch, ...
          fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;

  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    inputs = params.getBatch(params.imdb, batch) ;

    if params.prefetch
      if s == params.numSubBatches
        batchStart = t + (labindex-1) + params.batchSize ;
        batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
      params.getBatch(params.imdb, nextBatch) ;
    end

    if isa(net, 'dagnn.DagNN')
      if strcmp(mode, 'train')
        net.mode = 'normal' ;
        net.accumulateParamDers = (s ~= 1) ;
        net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
      else
        net.mode = 'test' ;
        net.eval(inputs) ;
      end
    else  % AutoNN
      if strcmp(mode, 'train')
        net.eval(inputs, 'normal', params.derOutputs, s ~= 1) ;
      else
        net.eval(inputs, 'test') ;
      end
    end
  end

  % Check if any vars contain NaN or Inf
  if ~isempty(params.onNanInf) && ~isa(net, 'dagnn.DagNN') && ...
    ~all(cellfun(@(x) ~isnumeric(x) || gather(all(isfinite(x(:)))), net.vars))
    switch params.onNanInf
    case 'error'
      error('The network contains NaN or Inf values.') ;
    case 'warning'
      warning('The network contains NaN or Inf values.') ;
    otherwise  % 'debug'/'break'
      net.displayVars() ;
      warning('The network contains NaN or Inf values.') ;
      keyboard ;
    end
  end

  % Accumulate gradient.
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    if isa(net, 'dagnn.DagNN')
      state = accumulateGradients(net, state, params, batchSize, parserv) ;
    else
      state = accumulateGradientsAutoNN(net, state, params, batchSize, parserv) ;
    end
  end

  % Get statistics.
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats.num = num ;
  stats.time = time ;
  stats = params.extractStatsFn(stats, net, batchSize / max(1, numGpus)) ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

  fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s: %.3f', f, stats.(f)) ;
  end
  fprintf('\n') ;

  if params.plotDiagnostics && ~isa(net, 'dagnn.DagNN') && mod(t-1, params.batchSize * 5) == 0
    net.plotDiagnostics(200) ;
  end
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveSolverState
  state.solverState = [] ;
else
  for i = 1:numel(state.solverState)
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gather(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gather, s, 'UniformOutput', false) ;
    end
  end
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;
if params.normalizedLoss  % assume loss is normalized by batch size
  batchSize = max(1,numGpus) * params.numSubBatches ;
end

for p=1:numel(net.params)

  if ~isempty(parserv)
    parDer = parserv.pullWithIndex(p) ;
  else
    parDer = net.params(p).der ;
  end

  switch net.params(p).trainMethod

    case 'average' % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      net.params(p).value = vl_taccum(...
          1 - thisLR, net.params(p).value, ...
          (thisLR/batchSize/net.params(p).fanout),  parDer) ;

    case 'gradient'
      thisDecay = params.weightDecay * net.params(p).weightDecay ;
      thisLR = params.learningRate * net.params(p).learningRate ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/batchSize, parDer, ...
                           thisDecay, net.params(p).value) ;

        if isempty(params.solver)
          % Default solver is the optimised SGD.
          % Update momentum.
          state.solverState{p} = vl_taccum(...
            params.momentum, state.solverState{p}, ...
            -1, parDer) ;

          % Nesterov update (aka one step ahead).
          if params.nesterovUpdate
            delta = params.momentum * state.solverState{p} - parDer ;
          else
            delta = state.solverState{p} ;
          end

          % Update parameters.
          net.params(p).value = vl_taccum(...
            1,  net.params(p).value, thisLR, delta) ;

        else
          % call solver function to update weights
          [net.params(p).value, state.solverState{p}] = ...
            params.solver(net.params(p).value, state.solverState{p}, ...
            parDer, params.solverOpts, thisLR) ;
        end
      end
    otherwise
      error('Unknown training method ''%s'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

% -------------------------------------------------------------------------
function state = accumulateGradientsAutoNN(net, state, params, batchSize, parserv)
% -------------------------------------------------------------------------

% ensure supported training methods are ordered as expected
assert(isequal(Param.trainMethods, {'gradient', 'average', 'none'})) ;

numGpus = numel(params.gpus) ;
if params.normalizedLoss  % assume loss is normalized by batch size
  batchSize = max(1,numGpus) * params.numSubBatches ;
end

paramVars = [net.params.var] ;
w = net.getValue(paramVars) ;
dw = net.getDer(paramVars) ;
if isscalar(paramVars), w = {w} ; dw = {dw} ; end

if ~params.plotDiagnostics
  % allow memory to be released, for parameters and their derivatives
  net.setValue([paramVars, paramVars + 1], cell(1, 2 * numel(paramVars))) ;
else
  % free only parameter memory, as we still need the gradients for plotting the diagnostics
  net.setValue(paramVars, cell(size(paramVars))) ;
end

for p=1:numel(net.params)
  if ~isempty(parserv)
    parDer = parserv.pullWithIndex(p) ;
  else
    parDer = dw{p} ;
  end

  switch net.params(p).trainMethod
    case 1
      thisDecay = params.weightDecay * net.params(p).weightDecay ;
      thisLR = params.learningRate * net.params(p).learningRate ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/batchSize, parDer, ...
                           thisDecay, w{p}) ;

        if isempty(params.solver)
          % Default solver is the optimised SGD.
          % Update momentum.
          state.solverState{p} = vl_taccum(...
            params.momentum, state.solverState{p}, ...
            -1, parDer) ;

          % Nesterov update (aka one step ahead).
          if params.nesterovUpdate
            delta = params.momentum * state.solverState{p} - parDer ;
          else
            delta = state.solverState{p} ;
          end

          % Update parameters.
          w{p} = vl_taccum(1, w{p}, thisLR, delta) ;

        else
          % call solver function to update weights
          [w{p}, state.solverState{p}] = ...
            params.solver(w{p}, state.solverState{p}, ...
            parDer, params.solverOpts, thisLR) ;
        end
      end

    case 2 % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      w{p} = vl_taccum(...
          1 - thisLR, w{p}, ...
          (thisLR/batchSize/net.params(p).fanout),  parDer) ;

    case 3  % none
    otherwise
      error('Unknown training method ''%i'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

if isscalar(paramVars), w = w{1} ; end
net.setValue(paramVars, w) ;

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net, sel, ~, ~)
% -------------------------------------------------------------------------
for i = 1:numel(sel)
  stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function stats = extractStatsAutoNN(stats, net, sel, batchSize, normalizedLoss)
% -------------------------------------------------------------------------
for i = 1:numel(sel)
  name = net.forward(sel(i)).name ;
  if ~isfield(stats, name)
    stats.(name) = 0 ;
  end
  newValue = gather(sum(net.vars{net.forward(sel(i)).outputVar(1)}(:))) ;
  if normalizedLoss  % Undo normalization for following code
    newValue = newValue * batchSize ;
  end
  % Update running average (same work as dagnn.Loss)
  stats.(name) = ((stats.num - batchSize) * stats.(name) + newValue) / stats.num ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
if isfield(net, 'layers')
%   net = dagnn.DagNN.loadobj(net) ;
  error('The last checkpoint does not contain a valid AutoNN model.') ;
else
  net = Net(net) ;
end
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end

end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename)
  clearMex() ;
  if numGpus == 1
    gpuDevice(opts.gpus)
  else
    spmd
      clearMex() ;
      gpuDevice(opts.gpus(labindex))
    end
  end
end


text:
%glinspace - Gamma-corrected linearly spaced vector.
%
%  USAGE
%
%    x = glinspace(d1,d2,gamma,n)
%
%    d1             start value
%    d2             stop value
%    gamma          gamma value
%    n              number of values
%
%  SEE ALSO
%
%    See also glinspace.
%

% Copyright (C) 2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

function x = glinspace(d1,d2,gamma,n)

% Check number of parameters
if nargin < 3,
  error('Incorrect number of parameters (type ''help glinspace'' for details).');
end

% Check parameters
if ~isdscalar(d1),
  error('Incorrect start value (type ''help glinspace'' for details).');
end
if ~isdscalar(d2),
  error('Incorrect stop value (type ''help glinspace'' for details).');
end
if ~isdscalar(gamma) || gamma <= 0,
  error('Incorrect start value (type ''help glinspace'' for details).');
end
if nargin < 3,
	n = 100;
end
if ~isiscalar(n),
  error('Incorrect number of values (type ''help glinspace'' for details).');
end

if d1 > d2
	gamma = 1/gamma;
end

x = linspace(0,1,n) .^ (1/gamma);
x = x *(d2-d1) + d1;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function out = issmooth(f)
%ISSMOOTH   True if a SINGFUN object F is smooth.
%   ISSMOOTH(F) returns TRUE if the SINGFUN objects F has zero EXPONENTS or if
%   the SMOOTHPART is zero. The test is FALSE otherwise.
%
% See also ISSING.

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

out = all(f.exponents == 0) || all(iszero(f.smoothPart));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function options = gpnddisimOptions

% GPNDDISIMOPTIONS Creates a set of default options for a GPNDDISIM model.
% FORMAT
% DESC returns a default options stucture for a GPNDDISIM model.
% RETURN options : the options structure.
%
% SEEALSO : gpnddisimCreate
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2008
%
% COPYRIGHT : Pei Gao, 2008
%
% COPYRIGHT :  Jaakko Peltonen, 2011

% GPNDDISIM

options=struct();
options.includeNoise=1;
options.addPriors=0;
options.optimiser='quasinew';
options.use_disimstartmean=1;
options.use_fixedrnavariance=0;
options.fix=[];

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function RV = fmri_vrestriction(nH, nC, AC, CC, nHTest)
%
% RV = fmri_vrestriction(nH, nC, AC, CC, nHTest)
%
% Creates a restriction vector of dimension 1 X nC*nH
%
% nH - total number of elements in HDIR
% nC - total number of stimulus conditions (including fixation)
% AC - list of active conditions
% CC - list of control conditions
% nHTest - list of HDIR components to test (default: all).
%
% Eg: fmri_vrestriction(10, 6, [1 3], [2 5], [3 6:9])
% Generates an RV for testing conditions (1+3)-(2+5) using
% components 3,6,7,8, and 9 in the HDIR.  Conditions 0 and 4
% are not tested. There are 10 components in the HDIR, and 5
% stimulus conditions.
%
%
% $Id: fmri_vrestriction.m,v 1.1 2004/03/12 07:31:56 sayres Exp $

% Check for the correct number of arguments %
if(nargin ~= 4 & nargin ~= 5)
  error('USAGE: fmri_vrestriction(nH, nC, AC, CC, <nHTest>');
  return;
end

%% check that the input parameters are correct %%
if(nH < 1)
  msg = 'nH must be greater than 0';
  qoe(msg); error(msg);
end
if(nC < 1)
  msg = 'Number of conditions must be greater than 0';
  qoe(msg); error(msg);
end
if(length(find(AC<0 | AC > nC-1)))
  msg = 'Invalid condition number in AC';
  qoe(msg); error(msg);
end
if(length(find(CC<0 | CC > nC-1)))
  msg = 'Invalid condition number in CC';
  qoe(msg); error(msg);
end

% Test that the same condition is not in both AC and CC %
C = [AC CC];
for n = 1:nC,
  if(length(find(C==n)) > 1)
    msg = 'Same condition is in both AC and CC';
    qoe(msg); error(msg);
  end
end

% Set defaults for nHTest, if needed, else test range %
if(nargin == 4) nHTest = [1:nH];
else
  if(length(find(nHTest<1 | nHTest > nH)))
    msg = 'Invalid nHTest Component';
    qoe(msg); error(msg);
  end
end

% Strip Condition 0 from AC (if there) %
iACnz = find(AC ~= 0);
if( ~ isempty(iACnz) )
  AC = AC(iACnz);
end
% Strip Condition 0 from CC (if there) %
iCCnz = find(CC ~= 0);
if( ~ isempty(iCCnz) )
  CC = CC(iCCnz);
end

%% -------- Generate the restriction matrix -------- %%
RV = zeros(nC-1,nH); % -1 excludes fixation
if( ~ isempty(iACnz) )
  RV(AC,nHTest) =  ones(length(AC),length(nHTest));
end
if( ~ isempty(iCCnz) )
  RV(CC,nHTest) = -ones(length(CC),length(nHTest));
end
RV = reshape(RV', [1 prod(size(RV))]);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function adj = meshAdjacencyMatrix(faces, varargin)
%MESHADJACENCYMATRIX Compute adjacency matrix of a mesh from set of faces.
%
%   ADJMAT = meshAdjacencyMatrix(FACES)
%   Returns a sparse NV-by-NV matrix (NV being the largest vertex index)
%   containing vertex adjacency of the mesh represented by FACES.
%   FACES is either a NF-by-3, a NF-by-4 index array, or a Nf-by-1 cell
%   array.
%
%   Example
%     [v f] = createCube;
%     adj = meshAdjacencyMatrix(f);
%
%   See also
%     meshes3d, triangulateFaces, smoothMesh

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2013-04-30,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2013 INRA - Cepia Software Platform.

% Ensures faces is a N-by-3 or N-by-4 array
if iscell(faces) || (isnumeric(faces) && size(faces, 2) > 4)
    faces = triangulateFaces(faces);
end

% forces faces to be floating point array, for sparse function
if ~isfloat(faces)
    faces = double(faces);
end
nv = max(faces(:));

% populate a sparse matrix
if size(faces, 2) == 3
    adj = sparse(...
        [faces(:,1); faces(:,1); faces(:,2); faces(:,2); faces(:,3); faces(:,3)], ...
        [faces(:,3); faces(:,2); faces(:,1); faces(:,3); faces(:,2); faces(:,1)], ...
        1.0, nv, nv);
elseif size(faces, 2) == 4
    adj = sparse(...
        [faces(:,1); faces(:,1); faces(:,2); faces(:,2); faces(:,3); faces(:,3); faces(:,4); faces(:,4)], ...
        [faces(:,4); faces(:,2); faces(:,1); faces(:,3); faces(:,2); faces(:,4); faces(:,3); faces(:,1)], ...
        1.0, nv, nv);
end

% remove double adjacencies
adj = min(adj, 1);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [c,bb] = hcs(riskscore,y,z,t,varargin)
%HCS Compute Harrell's C for survival model at given time
%
%  Description
%    [C, BB] = HCS(RISKSCORE,Y,Z,T,OPTIONS) Given a risk score vector
%    RISKSCORE, an observed time vector Y, a censoring indicator column
%    vector Z (0=event, 1=censored) and time T, returns
%    Harrell's C at time T and its estimated density using
%    Bayesian Bootstrap method. Large value of RISKSCORE should predict
%    early event.
%
%    The implemented estimator is called ExtAUC(t)_Count in the reference.
%
%    OPTIONS is optional parameter-value pair
%      rsubstream - number of a random stream to be used for
%                   simulating dirrand variables. This way same
%                   simulation can be obtained for different models.
%                   See doc RandStream for more information.
%
%  Reference
%    L. E. Chambless, C. P. Cummiskey, and G. Cui (2011). Several
%    methods to assess improvement in risk prediction models:
%    Extension to survival analysis. Statistics in Medicine
%    30(1):22-38.

% Copyright (C) 2012 Tomi Peltola, Ernesto Ulloa, Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

ip=inputParser;
ip.addRequired('riskscore',@(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) isreal(x) && all(isfinite(x(:))))
ip.addRequired('z', @(x) isreal(x) && all(isfinite(x(:))))
ip.addRequired('t', @(x) isreal(x) && isscalar(x) && ~isnan(x))
ip.addParamValue('rsubstream',0,@(x) isreal(x) && isscalar(x) && isfinite(x) && x>0)
ip.parse(riskscore,y,z,t,varargin{:})
rsubstream=ip.Results.rsubstream;

n=size(riskscore,1);

if nargout < 2
    % no bootstrapping
    numer = 0;
    denom = 0;
    for i = 1:n
        % i should have had event before or at t
        if z(i) == 1 || y(i) > t
            continue;
        end
        for j = 1:n
            % require y(i) < y(j) and check if the pair is concordant
            if y(i) >= y(j), continue, end;
            numer = numer + (riskscore(i) > riskscore(j));
            denom = denom + 1;
        end
    end
    c = numer / denom;
else
    % with bootstrapping
    n_replicates = 1000;

    % get substream
    if rsubstream > 0
        prevstream=setrandstream(0,'mrg32k3a');
        stream.Substream = rsubstream;
    end

    % bootstrap weights (first column is unit weights)
    w = [ones(n, 1) dirrand(n, n_replicates)];

    % set RandStream back to previous
    if rsubstream > 0
        setrandstream(prevstream);
    end

    % same algorithm as above, now with weights
    numer = zeros(1, size(w, 2));
    denom = zeros(1, size(w, 2));
    for i = 1:n
        if z(i) == 1 || y(i) > t
            continue;
        end
        for j = 1:n
            if y(i) >= y(j), continue, end;
            w_ = w(i, :) .* w(j, :);
            numer = numer + w_ * (riskscore(i) > riskscore(j));
            denom = denom + w_;
        end
    end
    bb = numer ./ denom;
    c = bb(1);  % first replicate is with unit weights
    bb(1) = []; % others are the bootstrap replicates
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function a = plu ( n, pivot )

%*****************************************************************************80
%
%% PLU returns the PLU matrix.
%
%  Discussion:
%
%    The PLU matrix has known P, L and U Gauss factors.
%
%  Example:
%
%    Input:
%
%      N = 5
%      PIVOT = ( 1, 3, 3, 5, 5 )
%
%    Output:
%
%      A:
%
%         11            12           13            14           15
%          1.375         9.75        43.25         44.75        46.25
%          2.75         25           26.25         27.5         28.75
%          0.34375       2.4375       7.71875      17.625       73.125
%          0.6875        4.875       15.4375       60           61.5625
%
%      P:
%
%          1             0            0             0            0
%          0             0            1             0            0
%          0             1            0             0            0
%          0             0            0             0            1
%          0             0            0             1            0
%
%      L:
%
%         1              0            0             0            0
%         0.25           1            0             0            0
%         0.125          0.375        1             0            0
%         0.0625         0.1875       0.3125        1            0
%         0.03125        0.09375      0.15625       0.21875      1
%
%      U:
%
%        11             12           13            14           15
%         0             22           23            24           25
%         0              0           33            34           35
%         0              0            0            44           45
%         0              0            0             0           55
%
%  Note:
%
%    The LINPACK routine DGEFA will factor the above A as:
%
%       11             12             13             14             15
%      -0.125          22             23             24             25
%      -0.25           -0.375         33             34             35
%      -0.03125        -0.09375       -0.15625       44             45
%      -0.0625         -0.1875        -0.3125        -0.21875       55
%
%    and the pivot information in the vector IPVT as:
%
%      ( 1, 3, 3, 5, 5 ).
%
%    The LAPACK routine DGETRF will factor the above A as:
%
%      11              12             13             14             15
%      0.25            22             23             24             25
%      0.125            0.375         33             34             35
%      0.0625           0.1875         0.3125        44             45
%      0.03125          0.09375        0.15625        0.21875       55
%
%   and the pivot information in the vector IPIV as:
%
%     ( 1, 3, 3, 5, 5 ).
%
%  Method:
%
%    The L factor will have unit diagonal, and subdiagonal entries
%    L(I,J) = ( 2 * J - 1 ) / 2^I, which should result in a unique
%    value for every entry.
%
%    The U factor of A will have entries
%    U(I,J) = 10 * I + J, which should result in "nice" entries as long
%    as N < 10.
%
%    The P factor can be deduced by applying the pivoting operations
%    specified by PIVOT in reverse order to the rows of the identity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer PIVOT(N), the list of pivot rows.  PIVOT(I)
%    must be a value between I and N, reflecting the choice of
%    pivot row on the I-th step.  For no pivoting, set PIVOT(I) = I.
%
%    Output, real A(N,N), the matrix.
%
  [ p, l, u ] = plu_plu ( n, pivot );

  a = p * l * u;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mutationChildren = mutationNcK(parents,options,GenomeLength,FitnessFcn,state,thisScore,thisPopulation)
% Oren Rosen
% The MathWorks
% 8/29/2007
%
% This custom mutation function is written to work on a population of
% vectors of zeros and ones with the same amount of ones in each vector.
% The mutated child is formed by randomly permuting the elements of the
% parent.
% Note: Performance-wise this hasn't worked out to be that efficient. A
% better implementation may swap only two of the elements.

mutationChildren = zeros(length(parents),GenomeLength);
numVars = length(thisPopulation(1,:));

for i=1:length(parents)
    child = thisPopulation(parents(i),:);
    mutationChildren(i,:) = child( randperm(numVars) );
end

author: Sable

repo: mcbench-benchmarks


text:
function output = callbonmin(model)

model = yalmip2nonlinearsolver(model);
options = [];
try
    options.bonmin = optiRemoveDefaults(model.options.bonmin,bonminset());
catch
    options.bonmin = model.options.bonmin;
end
options.ipopt = model.options.ipopt;
options.display = model.options.verbose;

if ~model.derivative_available
    disp('Derivate-free call to bonmin/ipopt not yet implemented')
    error('Derivate-free call to bonmin/ipopt not yet implemented')
end

if model.options.savedebug
    save bonmindebug model
end

Fupp = [ repmat(0,length(model.bnonlinineq)+length(model.K.q)*(model.K.q(1)>0),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(0,length(model.b),1);
    repmat(0,length(model.beq),1)];

Flow = [ repmat(-inf,length(model.bnonlinineq)+length(model.K.q)*(model.K.q(1)>0),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(-inf,length(model.b),1);
    repmat(0,length(model.beq),1)];

if isempty(Flow)
    Flow = [];
    Fupp = [];
end

% Since ipopt react strangely on lb>ub, we should bail if that is detected
% (ipopt creates an exception)
if ~isempty(model.lb)
    if any(model.lb>model.ub)
        problem = 1;
        solverinput = [];
        solveroutput = [];
        output = createoutput(model.lb*0,[],[],problem,'BONMIN',solverinput,solveroutput,0);
        return
    end
end

% These are needed to avoid recomputation due to ipopts double call to get
% f and df, and g and dg
global latest_x_f
global latest_x_g
global latest_df
global latest_f
global latest_G
global latest_g
global latest_xevaled
global latest_x_xevaled
latest_G= [];
latest_g = [];
latest_x_f = [];
latest_x_g = [];
latest_xevaled = [];
latest_x_xevaled = [];

funcs.objective = @(x)ipopt_callback_f(x,model);
funcs.gradient = @(x)ipopt_callback_df(x,model);
if ~isempty(Fupp)
    funcs.constraints = @(x)ipopt_callback_g(x,model);
    funcs.jacobian  = @(x)ipopt_callback_dg(x,model);
end

options.lb = model.lb(:)';
options.ub = model.ub(:)';
if ~isempty(Fupp)
    options.cl = Flow;
    options.cu = Fupp;
end

if ~isempty(Fupp)
    m = length(model.lb);
    allA=[model.Anonlinineq];
    if size(model.F_struc,1) > 0
        % These are SOCP cones
        top = 1;
        for i = 1:length(model.K.q)
            rows = model.F_struc(top:top + model.K.q(i)-1,2:end)
            allA = [allA;any(rows,1)];
        end
    end
    allA = [allA;model.Anonlineq];
    jacobianstructure = spalloc(size(allA,1),m,0);
    depends = allA | allA;
    for i = 1:size(depends,1)
        vars = find(depends(i,:));
        [ii,vars] = find(model.deppattern(vars,:));
        vars = unique(vars);
        s = size(jacobianstructure,1);
        for j = 1:length(vars)
            jacobianstructure(i,find(vars(j) == model.linearindicies)) = 1;
        end
    end
    allA=[model.A; model.Aeq];
    depends = allA | allA;
    jacobianstructure = [jacobianstructure;depends];

    Z = sparse(jacobianstructure);
    funcs.jacobianstructure = @() Z;
end

if ~model.options.usex0
    model.x0 = (options.lb+options.ub)/2;
    model.x0(isinf(options.ub)) = options.lb(isinf(options.ub))+1;
    model.x0(isinf(options.lb)) = options.ub(isinf(options.lb))-1;
    model.x0(isinf(model.x0)) = 0;
end

if ~isempty(model.binary_variables) | ~isempty(model.integer_variables)
    options.var_type = zeros(length(model.linearindicies),1);
    options.var_type(model.binary_variables) = -1;
    options.var_type(model.integer_variables) = 1;
end

showprogress('Calling BONMIN',model.options.showprogress);
solvertime = tic;
[xout,info] = bonmin(model.x0,funcs,options);
solvertime = toc(solvertime);

x = RecoverNonlinearSolverSolution(model,xout);

switch info.status
    case {0,2}
        problem = 0;
    case 1
        problem = 1;
    case {-1}
        problem = 3;
    case {3}
        problem = 15;
    otherwise
        problem = -1;
end

% Duals currently not supported
D_struc = [];

% Save all data sent to solver?
if model.options.savesolverinput
    solverinput.x0 = model.x0;
    solverinput.model = model;
    solverinput.options = options;
else
    solverinput = [];
end

% Save all data from the solver?
if model.options.savesolveroutput
    solveroutput.x = xout;
    solveroutput.info = info;
else
    solveroutput = [];
end

% Standard interface
output = createoutput(x,D_struc,[],problem,'BONMIN',solverinput,solveroutput,solvertime);

% Code supplied by Jonatan Currie
function opts = removeDefaults(opts,defs)
oFn = fieldnames(opts);
for i = 1:length(oFn)
    label = oFn{i};
    if(isfield(defs,label))
        if(ischar(opts.(label)))
            if(strcmpi(defs.(label),opts.(label)))
                opts = rmfield(opts,label);
            end
        else
            if(defs.(label) == opts.(label))
                opts = rmfield(opts,label);
            end
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function f = p12_f ( n, x )

%*****************************************************************************80
%
%% P12_F evaluates the objective function for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  f = 0.0;
  for i = 1 : 99
    arg = i / 100.0;
    r = abs ( ( - 50.0 * log ( arg ) )^( 2.0 / 3.0 ) + 25.0 - x(2) );

    t = exp ( - r^x(3) / x(1) ) - arg;

    f = f + t * t;

  end

  return
end


text:
function plot_frame(pt,R)
framesize = 1;
x_axis = R(:,1)*framesize;
y_axis = R(:,2)*framesize;
z_axis = R(:,3)*framesize;
plot3([pt(1) x_axis(1)+pt(1)],[pt(2) x_axis(2)+pt(2)], [pt(3) x_axis(3)+pt(3)], 'b')
%hold on
plot3([pt(1) y_axis(1)+pt(1)],[pt(2) y_axis(2)+pt(2)], [pt(3) y_axis(3)+pt(3)], 'g')
plot3([pt(1) z_axis(1)+pt(1)],[pt(2) z_axis(2)+pt(2)], [pt(3) z_axis(3)+pt(3)], 'r')
%hold off

author: yorgoon

repo: minimum-snap-geometric-control


text:
function option_num = p00_option_num ( problem )

%*****************************************************************************80
%
%% P00_OPTION_NUM returns the number of options available for a problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Output, integer OPTION_NUM, the number of options available for
%    this problem.  OPTION_NUM is always at least 1.
%
  if (  problem == 1 )
    option_num = p01_option_num ( );
  elseif ( problem == 2 )
    option_num = p02_option_num ( );
  elseif ( problem == 3 )
    option_num = p03_option_num ( );
  elseif ( problem == 4 )
    option_num = p04_option_num ( );
  elseif ( problem == 5 )
    option_num = p05_option_num ( );
  elseif ( problem == 6 )
    option_num = p06_option_num ( );
  elseif ( problem == 7 )
    option_num = p07_option_num ( );
  elseif ( problem == 8 )
    option_num = p08_option_num ( );
  elseif ( problem == 9 )
    option_num = p09_option_num ( );
  elseif ( problem == 10 )
    option_num = p10_option_num ( );
  elseif ( problem == 11 )
    option_num = p11_option_num ( );
  elseif ( problem == 12 )
    option_num = p12_option_num ( );
  elseif ( problem == 13 )
    option_num = p13_option_num ( );
  elseif ( problem == 14 )
    option_num = p14_option_num ( );
  elseif ( problem == 15 )
    option_num = p15_option_num ( );
  elseif ( problem == 16 )
    option_num = p16_option_num ( );
  elseif ( problem == 17 )
    option_num = p17_option_num ( );
  elseif ( problem == 18 )
    option_num = p18_option_num ( );
  elseif ( problem == 19 )
    option_num = p19_option_num ( );
  elseif ( problem == 20 )
    option_num = p20_option_num ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_OPTION_NUM - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized problem index = %d\n', problem );
    error ( 'P00_OPTION_NUM - Fatal error!' );
  end

  return
end


text:
function out = isLaue(s)
% check whether s is a Laue group

if ~isempty(s.LaueRef)

  out = s.LaueRef == s;

else

  if s.id > 0
    out = s.id == symmetry.pointGroups(s.id).LaueId;
  else
    out = any(s.rot(:) == rotation.inversion);
  end

  % store Laue group to speed up further checks
  if out
    s.LaueRef = s;
  else
    s.LaueRef = makeLaue(s);
  end

end

author: mtex-toolbox

repo: mtex

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


text:
function [ axis, angle ] = rotation_mat2axis_3d ( a )

%*****************************************************************************80
%
%% ROTATION_MAT2AXIS_3D converts a rotation from matrix to axis format in 3D.
%
%  Discussion:
%
%    The computation is based on the fact that a rotation matrix must
%    have an eigenvector corresponding to the eigenvalue of 1, hence:
%
%      ( A - I ) * v = 0.
%
%    The eigenvector V is the axis of rotation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jack Kuipers,
%    Quaternions and Rotation Sequences,
%    Princeton, 1998.
%
%  Parameters:
%
%    Input, real A(3,3), the rotation matrix.
%
%    Output, real AXIS(3), the axis vector which remains
%    unchanged by the rotation.
%
%    Output, real ANGLE, the angular measurement of the
%    rotation about the axis, in radians.
%
  dim_num = 3;
%
%  Compute the normalized axis of rotation.
%
  axis(1) = a(3,2) - a(2,3);
  axis(2) = a(1,3) - a(3,1);
  axis(3) = a(2,1) - a(1,2);

  axis_norm = sqrt ( sum ( axis(1:dim_num).^2 ) );

  if ( axis_norm == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ROTATION_MAT2AXIS_3D - Fatal error!\n' );
    fprintf ( 1, '  A is not a rotation matrix,\n' );
    fprintf ( 1, '  or there are multiple axes of rotation.\n' );
    error ( 'ROTATION_MAT2AXIS_3D - Fatal error!' );
  end

  axis(1:dim_num) = axis(1:dim_num) / axis_norm;
%
%  Find the angle.
%
  angle = r8_acos ( 0.5 * ( a(1,1) + a(2,2) + a(3,3) - 1.0 ) );

  return
end


text:
function [surface_original] = ft_surfacerealign(cfg, surface_original)

% FT_SURFACEREALIGN realigns surface
%
% FIDUCIAL - You can apply a rigid body realignment based on three fiducial
% locations. After realigning, the fiducials in the input surface
% (typically nose, left and right ear) are along the same axes as the
% fiducials in the template surface set.
%
%   cfg.method         = string representing the method for aligning the surface
%                        'fiducial'        realign using three fiducials
%                        (e.g. NAS, LPA and RPA)
%
% If you want to realign the surface using fiducials, the target and the
% objective have to contain the three fiducials which relate , e.g.
%   cfg.target.elecpos(1,:)    = [110 0 0]     % location of the nose
%   cfg.target.elecpos(2,:)    = [0  90 0]     % location of the left ear
%   cfg.target.elecpos(3,:)    = [0 -90 0]     % location of the right ear
%   cfg.target.label       = {'NAS', 'LPA', 'RPA'}
%   cfg.objective.elecpos(1,:) = [0 -110 0]      % location of the nose
%   cfg.objective.elecpos(2,:) = [90   0 0]      % location of the left ear
%   cfg.objective.elecpos(3,:) = [-90  0 0]      % location of the right ear
%   cfg.objective.label    = {'NAS', 'LPA', 'RPA'}
%
% See also FT_MESHREALIGN

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

% DEPRECATED by roboos on 22 August 2017
% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1830
% support for this functionality can be removed mid 2018
warning('FT_SURFACEREALIGN is deprecated, please use FT_MESHREALIGN instead.')

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision              = '$Id$';
ft_nargin                = nargin;
ft_nargout               = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    surface_original
ft_preamble provenance surface_original

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

usefiducial              = isfield(cfg, 'target') & isfield(cfg, 'objective') & strcmp(cfg.method,'fiducial');

if usefiducial
  cfg.elec                 = cfg.objective;
  %use ft_electroderealign for realign the surface
  surface_realigned           = ft_electroderealign(cfg);
  %use transformation obtained by ft_electroderealign
  transform                = [surface_original.pos, ones(length(surface_original.pos),1)]*surface_realigned.homogeneous';
  surface_original.pos        = transform(:,1:3);
  surface_original.cfg        = surface_realigned.cfg;
else
  ft_error('Cannot perform ft_surfacerealign. Please read help ft_surfacerealign, and check your cfg')
end


text:
classdef Scale < dagnn.ElementWise
  properties
    size
    hasBias = true
  end

  methods

    function outputs = forward(obj, inputs, params)
      args = horzcat(inputs, params) ;
      outputs{1} = bsxfun(@times, args{1}, args{2}) ;
      if obj.hasBias
        outputs{1} = bsxfun(@plus, outputs{1}, args{3}) ;
      end
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      args = horzcat(inputs, params) ;
      sz = [size(args{2}) 1 1 1 1] ;
      sz = sz(1:4) ;
      dargs{1} = bsxfun(@times, derOutputs{1}, args{2}) ;
      dargs{2} = derOutputs{1} .* args{1} ;
      for k = find(sz == 1)
        dargs{2} = sum(dargs{2}, k) ;
      end
      if obj.hasBias
        dargs{3} = derOutputs{1} ;
        for k = find(sz == 1)
          dargs{3} = sum(dargs{3}, k) ;
        end
      end
      derInputs = dargs(1:numel(inputs)) ;
      derParams = dargs(numel(inputs)+(1:numel(params))) ;
    end

    function params = initParams(obj,sc)
      params{1} = ones(obj.size,'single') * sc ;
      if obj.hasBias
        params{2} = zeros(obj.size,'single') ;
      end
    end

    function obj = Scale(varargin)
      obj.load(varargin) ;
    end
  end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%{
    Copyright (c) 2015, Philipp Krähenbühl
    All rights reserved.

    Redistribution and use in source and binary forms, with or without
    modification, are permitted provided that the following conditions are met:
        * Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.
        * Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.
        * Neither the name of the Stanford University nor the
        names of its contributors may be used to endorse or promote products
        derived from this software without specific prior written permission.

    THIS SOFTWARE IS PROVIDED BY Philipp Krähenbühl ''AS IS'' AND ANY
    EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
    WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
    DISCLAIMED. IN NO EVENT SHALL Philipp Krähenbühl BE LIABLE FOR ANY
    DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
    (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
	 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
	 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
	 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%}

% Usage example:
% I = imread('peppers.png');
% os = OverSegmentation( I );
% imagesc( os.s() ) % To visualize the over-segmentation

classdef LpoOverSegmentation
    properties (SetAccess = private, GetAccess = public)
        c_s
    end
    methods
        function obj = LpoOverSegmentation( I, NS, thick_bnd, thin_bnd )
            if nargin == 0
                error( 'Input image required' );
            end
            if nargin < 2
                NS = 1000;
            end
            if nargin == 3
                error( 'Cannot have a thick boundary map without a thinned out one!' );
            end
            if nargin < 3
                obj.c_s = lpo_mex( 'newImageOverSegmentation', I, NS );
            else
                obj.c_s = lpo_mex( 'newImageOverSegmentation', I, NS, thick_bnd, thin_bnd );
            end
        end
        function sobj = saveobj(this)
            sobj = this;
            sobj.c_s = lpo_mex( 'ImageOverSegmentation_serialize', this.c_s );
        end
        function r = s(this)
            r = lpo_mex( 'ImageOverSegmentation_s', this.c_s );
        end
        function r = boundaryMap(this)
            r = lpo_mex( 'ImageOverSegmentation_boundaryMap', this.c_s );
        end
        function r = maskToBox(this,mask)
            r = lpo_mex( 'ImageOverSegmentation_maskToBox', this.c_s, mask );
        end
        function delete(this)
            lpo_mex( 'freeImageOverSegmentation', this.c_s );
        end
    end
    methods (Static)
        function obj = loadobj(obj)
            data = obj.c_s;
            obj.c_s = lpo_mex( 'newImageOverSegmentationEmpty' );
            lpo_mex( 'ImageOverSegmentation_unserialize', obj.c_s, data )
        end
    end
end

author: Cloud-CV

repo: object-proposals


text:
classdef stressgrid < ZmapHGridFunction
    % STRESSGRID calculate stress grid for event that have Dip, DipDirection, and Rake
    properties
        calcmethod
    end

    properties(Constant)
        PlotTag         = 'stressgrid';
        ReturnDetails   = cell2table({ ...
            ...
            'S1Trend',  'S1Trend','';...
            'S1Plunge', 'S1Plunge','';...
            'S2Trend',  'S2Trend','';...
            'S2Plunge', 'S2Plunge','';...
            'S3Trend',  'S3Trend','';...
            'S3Plunge', 'S3Plunge','';...
            'Variance', 'Variance',''...
            }, 'VariableNames', {'Names','Descriptions','Units'});

        CalcFields = {...
            'S1Trend',  'S1Plunge',...
            'S2Trend',  'S2Plunge',...
            'S3Trend',  'S3Plunge',...
            'Variance'};

        ExtDir = fullfile(ZmapGlobal.Data.hodi, 'external');
        ParameterableProperties = [];

        %reference may not be the correct one
        References = "Michael AJ (1984) Determination of stress from slip data: faults and folds. J Geophys Res 89:11517–11526. doi:10.1029/JB089iB13p11517";
    end
    methods
        function obj = stressgrid(zap, varargin)
            % STRESSGRID calculate stress tensor using Michaels  code

            % NotImplemented: or Gephards code

            obj@ZmapHGridFunction(zap, 'S1Trend');
            report_this_filefun();

            obj.parseParameters(varargin);
            obj.StartProcess();

        end

        function InteractiveSetup(obj)
            %{
                dx = 0.1;
                dy = 0.1 ;
                ni = 50;
                ra = ZG.ra;
                Nmin = 0;
            %}

            zdlg = ZmapDialog();
            zdlg.AddPopup('calcmethod','Calc Method',{'Michaels Method'},1,...
                'Choose the only method that is available. (sorry, no other options)');

            obj.AddDialogOption(zdlg, 'EventSelector')
            zdlg.Create('Name', 'Stress grid options', 'WriteToObj',obj, 'OkFcn',@obj.doIt);
        end

        function results = Calculate(obj)
            % get the grid-size interactively and
            % calculate the b-value in the grid by sorting
            % the seismicity and selectiong the ni neighbors
            % to each grid point

            % get new grid if needed
            d = pwd;
            slick_cmd = fullfile(obj.ExtDir, append_system_specific_postfix('slick'))
            slfast_cmd = fullfile(obj.ExtDir, append_system_specific_postfix('slfast'))

            try
                cd(obj.ExtDir);
                itter=1;
                % calculate at all points
                obj.gridCalculations(@calculation_function);

                cd(d);

            catch ME

                cd(d);
                rethrow(ME)

            end

            if nargout
                results=obj.Result.values;
            end
            f=gcf;
            try
            obj.view_stressmap()
            catch ME
                warning(ME.message)
            end
            figure(f);

            function bvg = calculation_function(b)
                % returns: S1Trend S1Plunge S2Trend S2Plunge S3Trend S3Plunge Variance Radius b-value

                %estimate the completeness and b-value
                % Take the focal mechanism from actual catalog
                % tmpi-input: [dip direction (East of North), dip , rake (Kanamori)]
                %{
                tic
                itter=itter+1;
                inputFile = "dataz"+itter;
                outputFile = inputFile + ".oput";
                outputFileBoot = inputFile + ".slboot";

                % disp('writing raw inversion data to : ' + inputFile);
                % Create file for inversion
                fid = fopen(inputFile,'w');
                fprintf(fid,'%s \n', 'Inversion data');
                for i=1:b.Count
                    fprintf(fid,'%7.3f  %7.3f  %7.3f\n', b.DipDirection(i), b.Dip(i), b.Rake(i));
                end
                fclose(fid);

                % slick calculates the best solution for the stress tensor according to
                % Michael(1987): creates data2.oput
                % disp('slick')
                [mystatus, myresult] = system(slick_cmd + inputFile);

                % Get data from data2.oput

                [fBeta, fStdBeta, fTauFit, fAvgTau, fStdTau] = import_slickoput(outputFile);
                % Stress tensor inversion
                [mystatus, myresult] = system(slfast_cmd + inputFile); % --> writes to outputFileBoot

                sGetFile = fullfile(obj.ExtDir, outputFileBoot); %'data2.slboot'

                % Description of data2
                % Line 1: Variance S11 S12 S13 S22 S23 S33 => Variance and components of stress tensor (S = sigma)
                % Line 2: Phi S1t S1p S2t S2p S3t S3p   => Phi is relative size S3/S1, t=trend, p=plunge (other description)

                                % Result matrix
                % S1Trend, S1Plunge, S2Trend, S2Plunge, S3Trend, S3Plunge, Variance // ignored: Radius, bvalue

                d0=load(sGetFile, inputFile);
                bvg_old = [d0(2,2:7) d0(1,1)];

                if isfile(inputFile)
                    delete(inputFile);
                end
                if isfile(outputFile)
                    delete(outputFile);
                end
                if isfile(outputFileBoot)
                    delete(outputFileBoot);
                end

                toc
                disp('^-old   new-v');
                %}
                %tic

                %[fBeta2, fStdBeta2, fTauFit2, fAvgTau2, fStdTau2]=slick([b.DipDirection b.Dip b.Rake]);
                b = b.getAddon('MomentTensor');
                [bvg(7), ~, ~, bvg(1), bvg(2), bvg(3), bvg(4), bvg(5), bvg(6)]= slick(...
                    [b.DipDirection b.Dip b.Rake]);
                % toc

            end
        end

        function view_stressmap(obj)
            % view_stressmap
            % pulled into function due to its tight bonding from an m file
            % Author: S. Wiemer
            % updated: 19.05.2005, j.woessner@sed.ethz.ch
            % turned into function by Celso G Reyes 2017
            %FIXME this has never been updated

            values = obj.Result.values;

            ZG=ZmapGlobal.Data;
            SA = 1;
            SA2 = 3;

            % Matrix bvg contains:
            % bvg : [S1Trend S1Plunge S2Trend S2Plunge S3Trend S3Plunge Variance Radius b-value]
            % ste : [S1Plunge S1Trend+180 S2Plunge S2Trend+180 S3Plunge S3Trend+180 Variance];
            %ste = [bvg(:,2) bvg(:,1)+180  bvg(:,4) bvg(:,3)+180 bvg(:,6) bvg(:,5)+180 bvg(:,7) ];
            ste = [values.S1Plunge values.S1Trend+180  values.S2Plunge values.S2Trend+180 values.S3Plunge values.S3Trend+180 values.Variance ];
            sor = ste;
            % sor : [S1Plunge S1Trend+270 S2Plunge S2Trend+180 S3Plunge S3Trend+180 Variance];
            sor(:,SA*2) = sor(:,SA*2)+90;

            % Create matrices
            normlap2=NaN(height(values),1);

            valueMap=normlap2;%reshape(normlap2,length(yvect),length(xvect));
            %s11 = valueMap;

            % Create figure
            % fig=figure('visible','off');
            l_normal =  ste(:,1) > 52 &   ste(:,5) < 35 ;
            l_notnormal = l_normal < 1;
            ax=axes(figure);
            plq = quiver(ax,values.x(l_notnormal),values.y(l_notnormal),-cos(sor(l_notnormal,SA*2)*pi/180),sin(sor(l_notnormal,SA*2)*pi/180),0.6,'.');
            set(plq,'LineWidth',0.5,'Color','k')
            px = get(plq,'Xdata');
            py = get(plq,'Ydata');

            % fig=figure('visible','off')
            ax=axes(figure);
            plq_n = quiver(ax,values.x(l_normal),values.y(l_normal),-cos(sor(l_normal,SA2*2)*pi/180),sin(sor(l_normal,SA2*2)*pi/180),0.6,'.');
            set(plq_n,'LineWidth',0.5,'Color','r')

            drawnow
            px_n = get(plq_n,'Xdata');
            py_n = get(plq_n,'Ydata');
            %close

            fig=figure('Name','Faulting style map','pos',[100 100 860 600]);
            watchon;
            %whitebg(gcf);
            set(fig,'color','w');
            ax=axes('pos',[0.12 0.12 0.8 0.8]);
            set(ax,'NextPlot','add')
            n = 0;
            l0 = []; l1 = []; l2 = []; l3 = []; l4 = []; l5 = [];

            ste = sor(l_notnormal,:);
            for i = 1:3:length(px)-1
                n = n+1;
                j = jet;
                col = floor(ste(n,SA*2-1)/60*62)+1;
                if col > 64
                    col = 64;
                end
                pl = line(ax,px(i:i+1),py(i:i+1),'Linewidth',1,'Markersize',1);
                ax.NextPLot = 'add';

                dx = px(i)-px(i+1);
                dy = py(i) - py(i+1);
                pl2 = line(ax, px(i),py(i),'Marker','o','Markersize',0.1,'Linewidth',0.5);
                l0 = pl2;
                pl3 = line(ax,[px(i) px(i)+dx],[py(i) py(i)+dy],'Linewidth',1);

                % Select faulting style according to Zoback(1992)
                thecolor = [0 0 0];
                if ste(n,1) < 40  && ste(n,3) > 45  && ste(n,5) < 20
                    thecolor = [0.2 0.8 0.2];
                    l3 = pl;
                end
                if ste(n,1) < 20  && ste(n,3) > 45  && ste(n,5) < 40
                    thecolor =[0.2 0.8 0.2];
                    l3 = pl;
                end
                if ste(n,1) > 40  &&                 ste(n,1) < 52  && ste(n,5) < 20
                    thecolor = [1 0 1];
                    l2 = pl;
                end
                if ste(n,1) < 20  &&                 ste(n,5) > 40  &&  ste(n,5) <  52
                    thecolor = [0 1 1];
                    l4 = pl;
                end
                if ste(n,1) < 37  && ste(n,5) > 47
                    thecolor = [0 0 1];
                    l5 = pl;
                end
                set([pl pl2 pl3], 'color', thecolor);

            end
            %drawnow
            ste = sor(l_normal,:);
            n = 0;

            for i = 1:3:length(px_n)-1
                n = n+1;j = jet;
                col = floor(ste(n,SA*2-1)/60*62)+1;
                if col > 64 ; col = 64; end
                dx_n = px_n(i)-px_n(i+1);
                dy_n= py_n(i) - py_n(i+1);
                pl_n = plot(ax, px_n(i:i+1),py_n(i:i+1),'k','Linewidth',1,'Markersize',1,'color',[ 0 0 0  ] );
                set(ax,'NextPlot','add')
                dx = px_n(i)- px_n(i+1);
                dy = py_n(i) - py_n(i+1);
                pl2_n = plot(ax, px_n(i),py_n(i),'ko','Markersize',0.1,'Linewidth',0.5,'color',[0 0 0] );
                l0 = pl2;
                pl3_n = plot(ax, [px_n(i) px_n(i)+dx_n],[py_n(i) py_n(i)+dy_n],'k','Linewidth',1,'color',[0 0 0] );

                if ste(n,1) > 52  &&                 ste(n,5) < 35
                    set([pl_n pl3_n],'color','r');
                    set(pl2_n,'color','r');
                    l1 = pl_n;
                end
            end

            if isempty(l1)
                pl2 = plot(px,py,'kx','Linewidth',1,'color','r'); l1 = pl2; %set(l1,'visible','off');
            end
            if isempty(l2)
                pl2 = plot(px,py,'kx','Linewidth',1,'color','m'); l2 = pl2; %set(l2,'visible','off');
            end
            if isempty(l3)
                pl2 = plot(px,py,'kx','Linewidth',1,'color',[0.2 0.8 0.2] ); l3 = pl2; %set(l3,'visible','off');
            end
            if isempty(l4)
                pl2 = plot(px,py,'kx','Linewidth',1,'color','c' ); l4 = pl2;% set(l4,'visible','off');
            end
            if isempty(l5)
                pl2 = plot(px,py,'kx','Linewidth',1,'color','b' ); l5 = pl2; %set(l5,'visible','off');
            end
            if isempty(l0)
                l0 = plot(px,py,'kx','Linewidth',1,'color',[0 0 0 ] );  % set(l0,'visible','off');
            end

            try
                legend([l1 l2 l3 l4 l5 l0],'NF','NS','SS','TS','TF','U');
            catch
                disp('Legend could not be drawn')
            end

            % Figure settings
            set(ax,'NextPlot','add')
            axis(ax, 'equal')
            % zmap_update_displays();
            set(ax,'PlotBoxAspectRatio',[0.827 1 1])
            %axis(ax,[ s2_west s1_east s4_south s3_north])
            title(ax,sprintf('%s;  %g to %g', name, t0b, teb),'FontSize',ZmapGlobal.Data.fontsz.s,...
                'Color','k','FontWeight','normal');
            xlabel(ax,'Longitude ','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
            ylabel(ax,'Latitude ','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
            set(ax,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','normal',...
                'FontWeight','normal','LineWidth',1,...
                'Box','on','TickDir','out','Ticklength',[0.01 0.01])

            watchoff;

            % View the variance map
            valueMap = r;
            ZG.shading_style = 'interp';
            set(ax,'NextPlot','add')

            obj.add_grid_centers();

            for n=1:numel(obj.features)
                ft=obj.ZG.features(obj.features{n});
                copyobj(ft,ax);
            end

        end
    end
    methods(Static)

        function h=  AddMenuItem(parent, zapFcn, varargin)
            % create a menu item
            label = 'Map Stress Tensor';
            h = uimenu(parent, 'Label', label,...
                'MenuSelectedFcn', @(~,~)XYfun.stressgrid(zapFcn()),...
                varargin{:});
        end
    end % static methods
end
%{
function bvg = stressgrid
    % stressgrid create a grid (interactively), calculate stress tensor using Michaels or Gephards code
    %
    % Incoming data:
    % requires catalog contains
    %   * DipDirection (East of North)
    %   * Dip,
    %   * Rake (Kanamori Convention)
    %
    % original: Stefan Wiemer 1/95
    %
    % updated by: J. Woessner
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();
    fs=filesep;

    % get the grid parameter
    % initial values
    %

    dx = 0.1;
    dy = 0.1 ;
    ni = 50;
    ra = ZG.ra;
    Nmin = 0;

    zdlg = ZmapDialog();
    zdlg.AddPopup('calcmethod','Calc Method',{'Michaels Method'},1,...
        'Choose the only method that is available. (sorry, no other options)');

    zdlg.AddEventSelector('EventSelector', obj.EventSelector)
    [res, okpressed]=zdlg.Create('Name', 'Stress grid');
    if ~okpressed
        return
    end
    disp(res)
    Grid=ZG.gridopt;
    EventSelector=res.EventSelector;

    origdir=pwd;
    try
    my_calculate();
    catch ME
        cd(origdir)
        rethrow ME
    end

end
%}

author: CelsoReyes


text:
classdef Derivative < Algorithm

    properties (Access = public)
        order = 1;
        delta = 1;
        inPlaceComputation = true;
    end

    methods (Access = public)
        function obj = Derivative(order, delta)
            if nargin > 0
                obj.order = order;
                if nargin > 1
                    obj.delta = delta;
                end
            end
            obj.name = 'derivative';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kSignal;
        end

        function derivative = compute(obj, data)
            if size(data,1) == 1
                derivative = data;
            else
                if obj.order == 1
                    derivative = obj.computeFirstOrderDerivative(data);
                elseif obj.order == 2
                    derivative = obj.computeSecondOrderDerivative(data);
                else
                    derivative = [];
                end
            end
        end

        function derivative = computeFirstOrderDerivative(obj,data)
            n = length(data);
            derivative = zeros(n,1);
            for i = 2 : n
                derivative(i) = (data(i) - data(i-1)) / obj.delta;
            end
        end

        function derivative = computeSecondOrderDerivative(obj,data)
            n = length(data);
            derivative = zeros(n,1);

            derivative(1) = (data(2) - data(1)) / obj.delta;
            derivative(n) = (data(n) - data(n-1)) / obj.delta;

            deltaSquared = obj.delta * obj.delta;
            for i = 2 : n-1
                derivative(i) = (data(i-1) - data(i) + data(i+1)) / deltaSquared;
            end
        end

        function str = toString(obj)
            str = sprintf('%s_%d_%.2f',obj.name,obj.order,obj.delta);
        end

        function editableProperties = getEditableProperties(obj)
            property1 = Property('order',obj.order,1,2);
            property2 = Property('delta',obj.delta);
            editableProperties = [property1,property2];
        end

        function metrics = computeMetrics(obj,input)
            n = size(input,1);
            flops = 6 * obj.order * size(input,1);
            if obj.inPlaceComputation
                memory = 1;
            else
                memory = n * Constants.kSensorDataBytes;
            end
            outputSize = n * Constants.kSensorDataBytes;
            metrics = Metric(flops,memory,outputSize);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
%DEMO_WFBT  Auditory filterbanks built using filterbank tree structures
%
%   This demo shows two specific constructions of wavelet filterbank trees
%   using several M-band filterbanks with possibly different M (making the
%   thee non-homogenous) in order to approximate auditory frequency bands
%   and a musical scale respectively.  Both wavelet trees produce perfectly
%   reconstructable non-redundant representations.
%
%   The constructions are the following:
%
%      1) Auditory filterbank, approximately dividing the frequency band
%         into intervals reminisent of the bask scale.
%
%      2) Musical filterbank, approximately dividing the freq. band into
%         intervals reminicent of the well-tempered musical scale.
%
%   Shapes of the trees were taken from fig. 8 and fig. 9 from the refernece.
%   Sampling frequency of the test signal is 48kHz as used in the article.
%
%   .. figure::
%
%      Frequency responses of the auditory filterbank.
%
%      Both axes are in logarithmic scale.
%
%   .. figure::
%
%      TF plot of the test signal using the auditory filterbank.
%
%   .. figure::
%
%      Frequency responses of the musical filterbank.
%
%      Both axes are in logarithmic scale.
%
%   .. figure::
%
%      TF plot of the test signal using the musical filterbank.
%
%   References: cucla99

% Load a test signal and resample it to 48 kHz since such sampling
% rate is assumed in the reference.
f = dctresample(greasy,3*numel(greasy));
fs = 48000;

% Bark-like filterbank tree
% Creating tree depicted in Figure 8 in the reference.
w = wfbtinit({'cmband3',1});
w = wfbtput(1,0,'cmband6',w);
w = wfbtput(1,1,'cmband3',w);
w = wfbtput(2,0:1,'cmband5',w);
w = wfbtput(2,2:3,'cmband2',w);

% Convert to filterbank
[g,a] = wfbt2filterbank(w);

% Plot frequency responses
figure(1);
filterbankfreqz(g,a,2*2048,'plot','fs',fs,'dynrange',30,'posfreq','flog');

% Do the transform
[c,info] = wfbt(f,w);
disp('The reconstruction should be close to zero:')
norm(f-iwfbt(c,info))

figure(2);
plotwavelets(c,info,fs,'dynrange',60);

% Well-tempered musical scale filterbank tree
% Creating tree depicted in Figure 9 in the reference.
w2 = wfbtinit({'cmband4',1});
w2 = wfbtput(1,0:1,'cmband6',w2);
w2 = wfbtput(2,0:1,'cmband4',w2);
w2 = wfbtput(3,1:4,'cmband4',w2);

% Convert to filterbank
[g2,a2] = wfbt2filterbank(w2);
figure(3);
filterbankfreqz(g2,a2,2*2048,'plot','fs',fs,'dynrange',30,'posfreq','flog');

[c2,info2] = wfbt(f,w2);
disp('The reconstruction should be close to zero:')
norm(f-iwfbt(c2,info2))

figure(4);
plotwavelets(c2,info2,fs,'dynrange',60);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [pospeakind,negpeakind]=peakdetect(signal)

%	PEAKDETECT peak detection
%
%	[pospeakind,negpeakind]=peakdetect(signal)
%
%	The positive and negative polarity (concave down and up) peak index vectors are
%	generated from the signal vector and graphically displayed.  Positive and negative
%	polarity peaks occur at points of positive to negative and negative to positive
%	slope adjacency, respectively.  The typically rare contingencies of peaks
%	occurring at the lagging edges of constant intervals are supported.  Complex
%	signals are modified to the modulus of the elements.  If unspecified, the signal
%	vector is entered after the prompt from the keyboard.

%	Implemented using MATLAB 6.0.0
%
%	Examples:
%
%	ť [p,n]=peakdetect([-1 -1 0 1 0 1 0 -1 -1])
%
%	p =
%
%	     4     6
%
%	n =
%
%	     1     5     8
%
%	ť [p,n]=peakdetect(cos(2*pi*(0:999999)/500000))
%
%	p =
%
%	           1      500001     1000000
%
%	n =
%
%	      250001      750001
%
%	Copyright (c) 2001
%	Tom McMurray
%	mcmurray@teamcmi.com

%	if signal is not input, enter signal or return for empty outputs

if ~nargin
   signal=input('enter signal vector or return for empty outputs\n');
   if isempty(signal)
      pospeakind=[];
      negpeakind=[];
      return
   end
end
sizsig=size(signal);

%	while signal is unsupported, enter supported signal or return for empty outputs

while isempty(signal)|~isnumeric(signal)|~all(all(isfinite(signal)))...
      |length(sizsig)>2|min(sizsig)~=1
   signal=input(['signal is empty, nonnumeric, nonfinite, or nonvector:\nenter '...
         'finite vector or return for empty outputs\n']);
   if isempty(signal)
      pospeakind=[];
      negpeakind=[];
      return
   end
   sizsig=size(signal);
end

%	if signal is complex, modify to modulus of the elements

if ~isreal(signal)
   signal=abs(signal);
end

%	if signal is constant, return empty outputs

if ~any(signal-signal(1))
   pospeakind=[];
   negpeakind=[];
   disp('constant signal graph suppressed')
   return
end
sizsig1=sizsig(1);
lensig=sizsig1;

%	if signal is a row vector, modify to a column vector

if lensig==1
   signal=signal(:);
   lensig=sizsig(2);
end
lensig1=lensig-1;
lensig2=lensig1-1;

%	if signal length is 2, return max/min as positive/negative polarity peaks

if ~lensig2
   [sig,pospeakind]=max(signal);
   [sig,negpeakind]=min(signal);
   disp('2 element signal graph suppressed')
   return
end

%	generate difference signal

difsig=diff(signal);

%	generate vectors corresponding to positive slope indices

dsgt0=difsig>0;
dsgt00=dsgt0(1:lensig2);
dsgt01=dsgt0(2:lensig1);

%	generate vectors corresponding to negative slope indices

dslt0=difsig<0;
dslt00=dslt0(1:lensig2);
dslt01=dslt0(2:lensig1);

%	generate vectors corresponding to constant intervals

dseq0=difsig==0;
dseq01=dseq0(2:lensig1);
clear difsig

%	positive to negative slope adjacencies define positive polarity peaks

pospeakind=find(dsgt00&dslt01)+1;

%	negative to positive slope adjacencies define negative polarity peaks

negpeakind=find(dsgt01&dslt00)+1;

%	positive slope to constant interval adjacencies initiate positive polarity peaks

peakind=find(dsgt00&dseq01)+1;
lenpeakind=length(peakind);

%	determine positive polarity peak terminations

for k=1:lenpeakind
   peakindk=peakind(k);
   l=peakindk+1;

%	if end constant interval occurs, positive polarity peak exists

   if l==lensig
      pospeakind=[pospeakind;peakindk];

%	else l<lensig, determine next nonzero slope index

   else
      dseq0l=dseq0(l);
      while dseq0l&l<lensig1
         l=l+1;
         dseq0l=dseq0(l);
      end

%	if negative slope or end constant interval occurs, positive polarity peaks exist

      if dslt0(l)|dseq0l;
         pospeakind=[pospeakind;peakindk];
      end
   end
end

%	negative slope to constant interval adjacencies initiate negative polarity peaks

peakind=find(dslt00&dseq01)+1;
lenpeakind=length(peakind);
clear dseq01

%	determine negative polarity peak terminations

for k=1:lenpeakind
   peakindk=peakind(k);
   l=peakindk+1;

%	if end constant interval occurs, negative polarity peak exists

   if l==lensig
      negpeakind=[negpeakind;peakindk];

%	else l<lensig, determine next nonzero slope index

   else
      dseq0l=dseq0(l);
      while dseq0l&l<lensig1
         l=l+1;
         dseq0l=dseq0(l);
      end

%	if positive slope or end constant interval occurs, negative polarity peaks exist

      if dsgt0(l)|dseq0l;
         negpeakind=[negpeakind;peakindk];
      end
   end
end
clear dsgt0 peakind

%	if initial negative slope occurs, initial positive polarity peak exists

if dslt00(1)
   pospeakind=[1;pospeakind];

%	elseif initial positive slope occurs, initial negative polarity peak exists

elseif dsgt00(1)
   negpeakind=[1;negpeakind];

%	else initial constant interval occurs, determine next nonzero slope index

else
   k=2;
   dseq0k=dseq0(2);
   while dseq0k
      k=k+1;
      dseq0k=dseq0(k);
   end

%	if negative slope occurs, initial positive polarity peak exists

   if dslt0(k)
      pospeakind=[1;pospeakind];

%	else positive slope occurs, initial negative polarity peak exists

   else
      negpeakind=[1;negpeakind];
   end
end
clear dsgt00 dslt0 dslt00 dseq0

%	if final positive slope occurs, final positive polarity peak exists

if dsgt01(lensig2)
   pospeakind=[pospeakind;lensig];

%	elseif final negative slope occurs, final negative polarity peak exists

elseif dslt01(lensig2)
   negpeakind=[negpeakind;lensig];
end
clear dsgt01 dslt01

%	if peak indices are not ascending, order peak indices

if any(diff(pospeakind)<0)
   pospeakind=sort(pospeakind);
end
if any(diff(negpeakind)<0)
   negpeakind=sort(negpeakind);
end

%	if signal is a row vector, modify peak indices to row vectors

if sizsig1==1
   pospeakind=pospeakind.';
   negpeakind=negpeakind.';
end

%	plot signal peaks

plot(0:lensig1,signal,pospeakind-1,signal(pospeakind),'b^',negpeakind-1,...
   signal(negpeakind),'bv')
xlabel('Sample')
ylabel('Signal')
grid


text:
function this = MrDataNd_split_epoch(this)
%Tests split_epoch by creating artificial trials for fmri_short example
%
%   Y = MrUnitTest()
%   run(Y, 'MrDataNd_split_epoch');
%
% This is a method of class MrUnitTest.
%
% IN
%
% OUT
%
% EXAMPLE
%   MrDataNd_split_epoch
%
%   See also MrUnitTest

% Author:   Lars Kasper
% Created:  2019-03-25
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

doPlotRoi = 0;
doPlotTrialsActiveVoxel = true; % plot active voxel for all events

%% Load short fMRI example
pathExample = tapas_uniqc_get_path('examples');
%pathExample = fullfile(pathExample, 'nifti', 'rest');
%fileExample = fullfile(pathExample, 'fmri_short.nii');
pathExample = fullfile(pathExample, 'nifti', 'paradigm_visual');
fileExample = fullfile(pathExample, 'fmri.nii.gz');
fileBehav = fullfile(pathExample, 'behav.mat');
x = MrImage(fileExample);

TR = 3; % s

% update erroneous TR
x.dimInfo.set_dims('t', 'resolutions', TR, 'samplingWidths', TR, ...
    'firstSamplingPoint', 0);

%% Create mask for time series evaluation; central voxel in all slices

nSamplesXY = x.dimInfo.nSamples({'x', 'y'});

M = x.select('t',1);

%iActiveVoxelXY = round([nSamplesXY/2, nSamplesXY/2]);
idxActiveVoxel = [67 109 7]; % visual voxel
M.data(:) = 0;

% spherical ROI around active voxel
M.data(idxActiveVoxel(1), idxActiveVoxel(2), idxActiveVoxel(3)) = 1;
M = M.imdilate(strel('sphere', 5), '3D');

%% Split into epochs
% alternative [onsets, durations, names] = ...
%        get_multiple_conditions_visual('bevav.log')
% and use onsets{1} and durations{1} for the 10 blocks
load(fileBehav);
onsetTimes = relativeTimeBlockStartSeconds(1:4:end); % first blocks of same kind

newPeriStimulusOnsets = 10; % number of bins, if single number, duration of stimulation block: 14.925 s

y = x.split_epoch(onsetTimes, newPeriStimulusOnsets);

%% Plot all PST time courses from active voxel
if doPlotTrialsActiveVoxel
    stringTitle = sprintf('All trials from active voxel [%d %d %d]', ...
        idxActiveVoxel(1), idxActiveVoxel(2), idxActiveVoxel(3));
    figure('Name', stringTitle);
    ally = squeeze(y.data(idxActiveVoxel(1),idxActiveVoxel(2),idxActiveVoxel(3),:,:));
    plot(y.dimInfo.samplingPoints{4}, ally);
    hold all;
    meany = mean(ally,2);
    stdy = std(ally,0,2);
    plot(y.dimInfo.samplingPoints{4}, meany, 'g-', 'LineWidth',4)
    plot(y.dimInfo.samplingPoints{4}, meany-stdy,'k:', 'LineWidth',2)
    plot(y.dimInfo.samplingPoints{4}, meany+stdy,'k:', 'LineWidth',2)
    title(stringTitle);
end

%% evaluate time series: extract roi from both raw and epoched data
x.extract_rois(M);
x.compute_roi_stats();

y.extract_rois(M);
y.compute_roi_stats();

% plot with corresponding time vector
if doPlotRoi
    x.rois{1}.plot()
    y.rois{1}.plot()
end

%% Actual unit test

expSolution = 0;
actSolution.data = 0;

absTol = 10e-7;

%% verify equality of expected and actual solution
% import matlab.unittests to apply tolerances for objects
this.verifyEqual(actSolution.data, expSolution, 'absTol', absTol);

author: translationalneuromodeling

repo: tapas


text:
classdef KeplerianElementSetVariable < AbstractOrbitModelVariable
    %KeplerianElementSetVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) KeplerianElementSet = KeplerianElementSet.getDefaultElements();

        lb(1,6) double = 0;
        ub(1,6) double = 0;

        varSma(1,1) logical = false;
        varEcc(1,1) logical = false;
        varInc(1,1) logical = false;
        varRaan(1,1) logical = false;
        varArg(1,1) logical = false;
        varTru(1,1) logical = false;
    end

    methods
        function obj = KeplerianElementSetVariable(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.varSma)
                x(end+1) = obj.varObj.sma;
            end

            if(obj.varEcc)
                x(end+1) = obj.varObj.ecc;
            end

            if(obj.varInc)
                x(end+1) = obj.varObj.inc;
            end

            if(obj.varRaan)
                x(end+1) = obj.varObj.raan;
            end

            if(obj.varArg)
                x(end+1) = obj.varObj.arg;
            end

            if(obj.varTru)
                x(end+1) = obj.varObj.tru;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            useTf = obj.getUseTfForVariable();

            lb = obj.lb(useTf);
            ub = obj.ub(useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lb;
            ub = obj.lb;
        end

        function setBndsForVariable(obj, lb, ub)
            if(length(lb) == 6 && length(ub) == 6)
                obj.lb = lb;
                obj.ub = ub;
            else
                useTf = obj.getUseTfForVariable();

                obj.lb(useTf) = lb;
                obj.ub(useTf) = ub;
            end
        end

        function useTf = getUseTfForVariable(obj)
            useTf = [obj.varSma        obj.varEcc        obj.varInc        obj.varRaan        obj.varArg        obj.varTru];
        end

        function setUseTfForVariable(obj, useTf)
            obj.varSma  = useTf(1);
            obj.varEcc = useTf(2);
            obj.varInc  = useTf(3);
            obj.varRaan = useTf(4);
            obj.varArg  = useTf(5);
            obj.varTru = useTf(6);
        end

        function updateObjWithVarValue(obj, x)
            xInd = 1;

            if(obj.varSma)
                obj.varObj.sma = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varEcc)
                obj.varObj.ecc = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varInc)
                obj.varObj.inc = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varRaan)
                obj.varObj.raan = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varArg)
                obj.varObj.arg = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varTru)
                obj.varObj.tru = x(xInd);
                xInd = xInd + 1; %#ok<NASGU>
            end
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            if(evtNum > 0)
                subStr = sprintf('Event %i',evtNum);
            else
                subStr = varLocType;
            end

            nameStrs = {sprintf('%s SMA', subStr), ...
                        sprintf('%s Eccentricity', subStr), ...
                        sprintf('%s Inclination', subStr), ...
                        sprintf('%s RAAN', subStr), ...
                        sprintf('%s Argument of Periapsis', subStr), ...
                        sprintf('%s True Anomaly', subStr)};

            nameStrs = nameStrs(obj.getUseTfForVariable());
        end

        function varsStoredInRad = getVarsStoredInRad(obj)
            varsStoredInRad = [false false true true true true];
        end
    end
end


text:
function d = distanz(x,y,algtype)
% DISTANZ : calculates the distances between all vectors in x and y.
%
% usage:
%    d = distanz(x,y);
%
% inputs:
%    x      matrix with col vectors
%    y      matrix with col vectors (default == x)
%    algtype   the type of algorithm that is used (default==3)
%
% outputs:
%    d      distance matrix, not squared
%
% note:
%    part of the code is inspired by dist.m of the nntoolbox, other
%    part adapted from Francis Bach who took it from Roland
%    Bunschoten.
%
% sth * 19apr2002
% Adapted from create.m, originally written by
% (c) Stefan Harmeling, 2006

if exist('algtype')~=1 || isempty(algtype), algtype = 3; end
switch algtype
 case 1  % inspired by dist.m
  if exist('y')~=1 || isempty(y)
    % here comes code just for x
    [rx,cx] = size(x);
    d = zeros(cx,cx);
    nuller = zeros(cx,1);
    for c = 1:cx
      d(c,:) = sum((x-x(:,c+nuller)).^2,1);
    end
  else
    % here comes code for x and y
    [rx,cx] = size(x);
    [ry,cy] = size(y);
    if rx~=ry, error('x and y do not fit'), end
    d = zeros(cx,cy);
    if cx>cy
      nuller = zeros(cx,1);
      for c = 1:cy
	d(:,c) = sum((x-y(:,c+nuller)).^2,1)';
      end
    else
      nuller = zeros(cy,1);
      for c = 1:cx
	d(c,:) = sum((x(:,c+nuller)-y).^2,1);
      end
    end
  end

 case 2  % same as case 1, but with repmat instead of nuller
  if exist('y')~=1 || isempty(y)
    % here comes code just for x
    [rx,cx] = size(x);
    d = zeros(cx,cx);
    nuller = zeros(cx,1);
    for c = 1:cx
      d(c,:) = sum((x-repmat(x(:,c),[1 cx])).^2,1);
    end
  else
    % here comes code for x and y
    [rx,cx] = size(x);
    [ry,cy] = size(y);
    if rx~=ry, error('x and y do not fit'), end
    d = zeros(cx,cy);
    if cx>cy
      nuller = zeros(cx,1);
      for c = 1:cy
	d(:,c) = sum((x-repmat(y(:,c),[1 cx])).^2,1)';
      end
    else
      nuller = zeros(cy,1);
      for c = 1:cx
	d(c,:) = sum((repmat(x(:,c),[1 cy])-y).^2,1);
      end
    end
  end

 case 3  % inspired by Roland Bunschoten
  if exist('y')~=1 || isempty(y)
    % here comes code just for x
    cx = size(x,2);
    xx = sum(x.*x,1); xz = x'*x;
    d = abs(repmat(xx',[1 cx]) - 2*xz + repmat(xx,[cx 1]));
  else
    % here comes code for x and y
    [rx,cx] = size(x);
    [ry,cy] = size(y);
    if rx~=ry, error('x and y do not fit'), end
    xx = sum(x.*x,1); yy = sum(y.*y,1);  xy = x'*y;
    d = abs(repmat(xx',[1 cy]) + repmat(yy,[cx 1]) - 2*xy);
  end
end

d = sqrt(d);

repo: gspbox


text:
function out = spm_dartel_warp(job)
% Register images to template data.
% format spm_dartel_warp(job)
%
% The outputs are flow fields.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_dartel_warp.m 4064 2010-09-03 12:57:10Z john $

code = 2;
st = job.settings;
n1 = numel(job.images);
n2 = numel(job.images{1});
NF = struct('NI',[],'vn',[1 1]);
NF(n1,n2) = struct('NI',[],'vn',[1 1]);

for i=1:n1,
    if numel(job.images{i}) ~= n2,
        error('Incompatible number of images');
    end;
    for j=1:n2,
        [pth,nam,ext,num] = spm_fileparts(job.images{i}{j});
        NF(i,j).NI        = nifti(fullfile(pth,[nam ext]));
        num               = [str2num(num) 1 1];
        NF(i,j).vn        = num(1:2);
    end;
end;

dm = [size(NF(1,1).NI.dat) 1];
dm = dm(1:3);

spm_progress_bar('Init',n2,'Registering');
numits = 0;
for it=1:numel(st.param),
    numits = numits + st.param(it).its;
end

for i=1:n2,

    f = zeros([dm,n1],'single');
    for j=1:n1,
        vn = NF(j,i).vn;
        f(:,:,:,j) = single(NF(j,i).NI.dat(:,:,:,vn(1),vn(2)));
    end
    f(~isfinite(f)) = 0;

    [pth,nam,ext] = fileparts(NF(1,i).NI.dat.fname);
    fprintf('*** %s ***\n', nam);

    NU = NF(1,i).NI;
    [pth,nam,ext]   = fileparts(NU.dat.fname);
    NU.dat.fname = fullfile(pth,['u_' nam '.nii']);
    NU.dat.dim   = [dm 1 3];
    NU.dat.dtype = 'float32-le';
    NU.dat.scl_slope = 1;
    NU.dat.scl_inter = 0;
    NU.descrip = 'Flow Field';
    if exist(NU.dat.fname,'file'),
        fprintf('Continuing registration from pre-existing parameters (%s)\n', NU.dat.fname);

        u = NU.dat(:,:,:,1,:);
        u = single(squeeze(u));
    else
        u = zeros([dm,3],'single');
    end;

    it0 = 0;
    for it=1:numel(st.param),
        param = st.param(it);
        prm   = [st.rform, param.rparam, st.optim.lmreg, ...
                 st.optim.cyc, st.optim.its, param.K, code];

        % Load the template
        NG = nifti(strvcat(param.template{:}));

        g  = squeeze(single(NG.dat(:,:,:,:,:)));
        if ~all(size(g)==size(f)),
            error('Incompatible dimensions between images and template');
        end
        for j=1:param.its,
            it0 = it0 + 1;
            [u,ll] = dartel3(u,f,g,prm);
            fprintf('%d \t%g\t%g\t%g\t%g\n',...
                it0,ll(1),ll(2),ll(1)+ll(2),ll(3));

            spm_progress_bar('Set',i-1 + it0/numits);
        end
    end
    create(NU);
    NU.dat(:,:,:,1,:) = reshape(u,[dm 1 3]);
end;
spm_progress_bar('Clear');

n2 = numel(job.images{1});
out.files = cell(n2,1);
for j=1:n2,
    [pth,nam,ext,num] = spm_fileparts(job.images{1}{j});
    fname             = fullfile(pth,['u_' nam '.nii']);
    out.files{j}      = fname;
end;


text:
function [ result, seed ] = sphere01_triangle_quad_00 ( n, v1, v2, v3, f, ...
  seed )

%*****************************************************************************80
%
%% SPHERE01_TRIANGLE_QUAD_00: quadrature over a triangle on the unit sphere.
%
%  Discussion:
%
%    This is a Monte Carlo approach.
%
%    The integral is approximated by averaging the values at N random points,
%    multiplied by the area.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sample points.
%
%    Input, real V1(3), V2(3), V3(3), the XYZ coordinates of
%    the vertices of the triangle.
%
%    Input, function v = f ( x ), evaluates the integrand at X.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real RESULT, the approximate integral.
%
  area = sphere01_triangle_vertices_to_area ( v1, v2, v3  );

  [ vc, seed ] = sphere01_triangle_sample ( n, v1, v2, v3, seed );

  quad = 0.0;
  for j = 1 : n
    quad = quad + f ( vc(1:3,j) );
  end

  result = quad * area / n;

  return
end


text:
function value = c8_i ( )

%*****************************************************************************80
%
%% C8_I returns the value of the imaginary unit, i as a C8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, complex VALUE, the value of complex i.
%
  value = i;

  return
end


text:
function c = prod(a,dim)
%PROD         Implements  prod(a,dim)  for gradients
%
%   c = prod(a,dim)
%
% functionality as Matlab function prod for matrices, parameter dim optional
%

% written  10/16/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  [m n] = size(a);
  if nargin==1,
    if m==1
      dim=2;
    else
      dim=1;
    end
  end

  if dim==1
    c = ones(1,n);
    for i=1:m
%VVVV c = c .* a(i,:);
      s.type = '()'; s.subs = {i,':'}; c = c .* subsref(a,s);
%AAAA Matlab V5.2 bug fix
    end
  else
    c = ones(m,1);
    for i=1:n
%VVVV c = c .* a(:,i);
      s.type = '()'; s.subs = {':',i}; c = c .* subsref(a,s);
%AAAA Matlab V5.2 bug fix
    end
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function infoz(var1)
    % make an info window for various routines
    %
    %   A.Allmann

    global slapinfo p1

    report_this_filefun();

    if var1==1               %first call in P-Value window

        % This is the info window text
        %
        ttlStr='Help Window                                 ';

        hlpStr1=...
            ['                                                             '
            'No help available on this topic. Please                      '
            'consult the manual or contact stefan@giseis.alaska.edu       '
            '                                                             '];

    elseif var1 == 2                 % second help message

        ttlStr='Stress Tensor Inversion Resuls              ';

        hlpStr1=...
            ['                                                             '
            'This is a polar projection of the Stress Tensor Inversion    '
            'results, using Gepards and Forsyths (1984) technique.        '
            'Show in black are the best fitting stress directions,        '
            'The red. yellow and blue symbols indicate the 95% confidence '
            'region. The histogrom in the corner shows the distribution   '
            'of R values in this data-set.                                '
            '                                                             '];

    end

    % display the message
    zmaphelp(ttlStr,hlpStr1)
end