text:
function [net, stats] = cnn_train_binary(net, imdb, getBatch, varargin)
%CNN_TRAIN  An example implementation of SGD for training CNNs
%    CNN_TRAIN() is an example learner implementing stochastic
%    gradient descent with momentum to train a CNN. It can be used
%    with different datasets and tasks by providing a suitable
%    getBatch function.
%
%    The function automatically restarts after each training epoch by
%    checkpointing.
%
%    The function supports training on CPU or on one or more GPUs
%    (specify the list of GPU IDs in the `gpus` option).

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

opts.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.momentum = 0.9 ;
opts.randomSeed = 0 ;
opts.memoryMapFile = fullfile(tempdir, 'matconvnet.bin') ;
opts.profile = false ;

opts.conserveMemory = true ;
opts.backPropDepth = +inf ;
opts.sync = false ;
opts.cudnn = true ;
opts.errorFunction = 'multiclass' ;
opts.errorLabels = {} ;
opts.plotDiagnostics = false ;
opts.plotStatistics = true;
opts = vl_argparse(opts, varargin) ;
opts.numGpus = numel(opts.gpus);
opts.bwn = true;
opts.xnor = false;

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 isnan(opts.train), opts.train = [] ; end
if isnan(opts.val), opts.val = [] ; end

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

net = vl_simplenn_tidy(net); % fill in some eventually missing values
net.layers{end-1}.precious = 1; % do not remove predictions, used for error
vl_simplenn_display(net, 'batchSize', opts.batchSize) ;

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  for i=1:numel(net.layers)
    if isfield(net.layers{i}, 'weights')
      J = numel(net.layers{i}.weights) ;
      if ~isfield(net.layers{i}, 'learningRate')
        net.layers{i}.learningRate = ones(1, J, 'single') ;
      end
      if ~isfield(net.layers{i}, 'weightDecay')
        net.layers{i}.weightDecay = ones(1, J, 'single') ;
      end
    end
  end
end

% setup error calculation function
hasError = true ;
if isstr(opts.errorFunction)
  switch opts.errorFunction
    case 'none'
      opts.errorFunction = @error_none ;
      hasError = false ;
    case 'multiclass'
      opts.errorFunction = @error_multiclass ;
      if isempty(opts.errorLabels), opts.errorLabels = {'top1err', 'top5err'} ; end
    case 'binary'
      opts.errorFunction = @error_binary ;
      if isempty(opts.errorLabels), opts.errorLabels = {'binerr'} ; end
    otherwise
      error('Unknown error function ''%s''.', opts.errorFunction) ;
  end
end

state.getBatch = getBatch ;
stats = [] ;

% -------------------------------------------------------------------------
%                                                        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, stats] = loadState(modelPath(start)) ;
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.

  state.epoch = epoch ;
  state.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  state.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  state.val = opts.val(randperm(numel(opts.val))) ;
  state.imdb = imdb ;

  if numel(opts.gpus) <= 1
    [net,stats.train(epoch),prof] = process_epoch(net, state, opts, 'train') ;
    [~,stats.val(epoch)] = process_epoch(net, state, opts, 'val') ;
    if opts.profile
      profview(0,prof) ;
      keyboard ;
    end
  else
    spmd(numGpus)
      [net_, stats_.train, prof_] = process_epoch(net, state, opts, 'train') ;
      [~, stats_.val] = process_epoch(net_, state, opts, 'val') ;
      if labindex == 1, savedNet_ = net_ ; end
    end
    net = savedNet_{1} ;
    stats__ = accumulateStats(stats_) ;
    stats.train(epoch) = stats__.train ;
    stats.val(epoch) = stats__.val ;
    if opts.profile
      mpiprofile('viewer', [prof_{:,1}]) ;
      keyboard ;
    end
    clear net_ stats_ stats__ savedNet_ ;
  end

  % save
  if ~evaluateMode
    saveState(modelPath(epoch), net, stats) ;
  end

  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

% -------------------------------------------------------------------------
function err = error_multiclass(opts, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
[~,predictions] = sort(predictions, 3, 'descend') ;

% be resilient to badly formatted labels
if numel(labels) == size(predictions, 4)
  labels = reshape(labels,1,1,1,[]) ;
end

% skip null labels
mass = single(labels(:,:,1,:) > 0) ;
if size(labels,3) == 2
  % if there is a second channel in labels, used it as weights
  mass = mass .* labels(:,:,2,:) ;
  labels(:,:,2,:) = [] ;
end

m = min(5, size(predictions,3)) ;

error = ~bsxfun(@eq, predictions, labels) ;
err(1,1) = sum(sum(sum(mass .* error(:,:,1,:)))) ;
err(2,1) = sum(sum(sum(mass .* min(error(:,:,1:m,:),[],3)))) ;

% -------------------------------------------------------------------------
function err = error_binary(opts, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
error = bsxfun(@times, predictions, labels) < 0 ;
err = sum(error(:)) ;

% -------------------------------------------------------------------------
function err = error_none(opts, labels, res)
% -------------------------------------------------------------------------
err = zeros(0,1) ;

% -------------------------------------------------------------------------
function  [net_cpu,stats,prof] = process_epoch(net, state, opts, mode)
% -------------------------------------------------------------------------

% initialize empty momentum
if strcmp(mode,'train')
  state.momentum = {} ;
  for i = 1:numel(net.layers)
    if isfield(net.layers{i}, 'weights')
      for j = 1:numel(net.layers{i}.weights)
        state.layers{i}.momentum{j} = 0 ;
      end
    end
  end
end

% move CNN  to GPU as needed
numGpus = numel(opts.gpus) ;
if numGpus >= 1
  net = vl_simplenn_move(net, 'gpu') ;
end
if numGpus > 1
  mmap = map_gradients(opts.memoryMapFile, net, numGpus) ;
else
  mmap = [] ;
end

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

subset = state.(mode) ;
num = 0 ;
stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;
adjustTime = 0 ;
res = [] ;
error = [] ;

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

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

    [im, labels] = state.getBatch(state.imdb, batch) ;

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

    if numGpus >= 1
      im = gpuArray(im) ;
    end

    if strcmp(mode, 'train')
      dzdy = 1 ;
      evalMode = 'normal' ;
    else
      dzdy = [] ;
      evalMode = 'test' ;
    end
    net.layers{end}.class = labels ;
    if opts.bwn == true
        net = bwn_net(net);
    end
    if opts.xnor == true
        net = bwn_net(net);
    end
    % forward and backward
    res = vl_simplenn(net, im, dzdy, res, ...
                      'accumulate', s ~= 1, ...
                      'mode', evalMode, ...
                      'conserveMemory', opts.conserveMemory, ...
                      'backPropDepth', opts.backPropDepth, ...
                      'sync', opts.sync, ...
                      'cudnn', opts.cudnn) ;
    % bwn net or xnor net
    if opts.bwn == true
        [bnet,net] = bwn_net_back(net);
    end
    if opts.xnor == true
        [bnet,net] = bwn_net_back(net);
    end
    % accumulate errors

    error = sum([error, [...
      sum(double(gather(res(end).x))) ;
      reshape(opts.errorFunction(opts, labels, res),[],1) ; ]],2) ;
  end

  % accumulate gradient
  if strcmp(mode, 'train')
    if ~isempty(mmap)
      write_gradients(mmap, net) ;
      labBarrier() ;
    end
    [state, net] = accumulate_gradients(state, net, res, opts, batchSize, mmap) ;
  end

  % get statistics
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats = extractStats(net, opts, error / num) ;
  stats.num = num ;
  stats.time = time ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == opts.batchSize + 1
    % compensate for the first iteration, which is an outlier
    adjustTime = 2*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:', f) ;
    fprintf(' %.3f', stats.(f)) ;
  end
  fprintf('\n') ;

  % collect diagnostic statistics
  if strcmp(mode, 'train') && opts.plotDiagnostics
    switchfigure(2) ; clf ;
    diagn = [res.stats] ;
    diagnvar = horzcat(diagn.variation) ;
    barh(diagnvar) ;
    set(gca,'TickLabelInterpreter', 'none', ...
      'YTick', 1:numel(diagnvar), ...
      'YTickLabel',horzcat(diagn.label), ...
      'YDir', 'reverse', ...
      'XScale', 'log', ...
      'XLim', [1e-5 1]) ;
    drawnow ;
  end
end

if ~isempty(mmap)
  unmap_gradients(mmap) ;
end

if opts.profile
  if numGpus <= 1
    prof = profile('info') ;
    profile off ;
  else
    prof = mpiprofile('info');
    mpiprofile off ;
  end
else
  prof = [] ;
end

net_cpu = vl_simplenn_move(net, 'cpu') ;

% -------------------------------------------------------------------------
function [bnet,net] = bwn_net_back(net)
% -------------------------------------------------------------------------
n_layers = numel(net.layers);
bnet = net;
for ii = 1:n_layers
    if strcmp(net.layers{ii}.type,'conv')
        net.layers{ii}.weights{1} = net.layers{ii}.Bs{1};
        net.layers{ii}.weights{2} = net.layers{ii}.Bs{2};
    end
end

% -------------------------------------------------------------------------
function [net] = bwn_net(net)
% -------------------------------------------------------------------------
n_layers = numel(net.layers);
for ii = 1:n_layers
    if strcmp(net.layers{ii}.type,'conv')
        for jj = 1:numel(net.layers{ii}.alpha{1})
            temp_W = net.layers{ii}.weights{1};
            temp_b = net.layers{ii}.weights{2};
            temp_BW = ones(size(temp_W),'single');
            temp_Bb = ones(size(temp_b),'single');

            temp_BW(net.layers{ii}.weights{1} < 0) = -1;
            temp_Bb(net.layers{ii}.weights{2} < 0) = -1;

            net.layers{ii}.alpha{1,jj} = sum(temp_W(:))/numel(net.layers{ii}.weights{1});
            net.layers{ii}.alpha{2,jj} = sum(temp_b(:))/numel(net.layers{ii}.weights{1});

            net.layers{ii}.Bs{1} = temp_W;
            net.layers{ii}.Bs{2} = temp_b;

            net.layers{ii}.weights{1} = net.layers{ii}.alpha{1,jj}*temp_BW;
            net.layers{ii}.weights{2} = net.layers{ii}.alpha{2,jj}*temp_Bb;
        end
    end
end

% -------------------------------------------------------------------------
function [state, net] = accumulate_gradients(state, net, res, opts, batchSize, mmap)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for l=numel(net.layers):-1:1
  for j=1:numel(res(l).dzdw)

    % accumualte gradients from multiple labs (GPUs) if needed
    if numGpus > 1
      tag = sprintf('l%d_%d',l,j) ;
      for g = otherGpus
        tmp = gpuArray(mmap.Data(g).(tag)) ;
        res(l).dzdw{j} = res(l).dzdw{j} + tmp ;
      end
    end

    if j == 3 && strcmp(net.layers{l}.type, 'bnorm')
      % special case for learning bnorm moments
      thisLR = net.layers{l}.learningRate(j) ;
      net.layers{l}.weights{j} = ...
        (1 - thisLR) * net.layers{l}.weights{j} + ...
        (thisLR/batchSize) * res(l).dzdw{j} ;
    else
      % standard gradient training
      thisDecay = opts.weightDecay * net.layers{l}.weightDecay(j) ;
      thisLR = state.learningRate * net.layers{l}.learningRate(j) ;
      state.layers{l}.momentum{j} = opts.momentum * state.layers{l}.momentum{j} ...
        - thisDecay * net.layers{l}.weights{j} ...
        - (1 / batchSize) * res(l).dzdw{j} ;
      net.layers{l}.weights{j} = net.layers{l}.weights{j} + ...
        thisLR * state.layers{l}.momentum{j} ;
    end

    % if requested, collect some useful stats for debugging
    if opts.plotDiagnostics
      variation = [] ;
      label = '' ;
      switch net.layers{l}.type
        case {'conv','convt'}
          variation = thisLR * mean(abs(state.layers{l}.momentum{j}(:))) ;
          if j == 1 % fiters
            base = mean(abs(net.layers{l}.weights{j}(:))) ;
            label = 'filters' ;
          else % biases
            base = mean(abs(res(l+1).x(:))) ;
            label = 'biases' ;
          end
          variation = variation / base ;
          label = sprintf('%s_%s', net.layers{l}.name, label) ;
      end
      res(l).stats.variation(j) = variation ;
      res(l).stats.label{j} = label ;
    end
  end
end

% -------------------------------------------------------------------------
function mmap = map_gradients(fname, net, numGpus)
% -------------------------------------------------------------------------
format = {} ;
for i=1:numel(net.layers)
  for j=1:numel(net.layers(i).params)
    par = net.layers(i).params{j} ;
    format(end+1,1:3) = {'single', size(par), sprintf('l%d_%d',i,j)} ;
  end
end
format(end+1,1:3) = {'double', [3 1], 'errors'} ;
if ~exist(fname) && (labindex == 1)
  f = fopen(fname,'wb') ;
  for g=1:numGpus
    for i=1:size(format,1)
      fwrite(f,zeros(format{i,2},format{i,1}),format{i,1}) ;
    end
  end
  fclose(f) ;
end
labBarrier() ;
mmap = memmapfile(fname, ...
                  'Format', format, ...
                  'Repeat', numGpus, ...
                  'Writable', true) ;

% -------------------------------------------------------------------------
function write_gradients(mmap, net, res)
% -------------------------------------------------------------------------
for i=1:numel(net.layers)
  for j=1:numel(res(i).dzdw)
    mmap.Data(labindex).(sprintf('l%d_%d',i,j)) = gather(res(i).dzdw{j}) ;
  end
end

% -------------------------------------------------------------------------
function unmap_gradients(mmap)
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
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(net, opts, errors)
% -------------------------------------------------------------------------
stats.objective = errors(1) ;
for i = 1:numel(opts.errorLabels)
  stats.(opts.errorLabels{i}) = errors(i+1) ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net, stats)
% -------------------------------------------------------------------------
save(fileName, 'net', 'stats') ;

% -------------------------------------------------------------------------
function [net, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'stats') ;
net = vl_simplenn_tidy(net) ;

% -------------------------------------------------------------------------
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 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
  if exist(opts.memoryMapFile)
    delete(opts.memoryMapFile) ;
  end
end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename)
  if numGpus == 1
    gpuDevice(opts.gpus)
  else
    spmd, gpuDevice(opts.gpus(labindex)), end
  end
end

author: jiangqy


text:
function varargout = mesh(varargin)
%MESH   Mesh surface of a DISKFUN.
%   MESH is not supported for DISKFUN objects, and so returns an error.
%
% See also DISKFUN/SURF.

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

[varargout{1:nargout}] = mesh@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef LaunchVehicleExtremaTypeEnum < matlab.mixin.SetGet
    %LaunchVehicleExtremaTypeEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        Maximum('Maximum');
        Minimum('Minimum');
    end

    properties
        nameStr char = '';
    end

    methods
        function obj = LaunchVehicleExtremaTypeEnum(nameStr)
            obj.nameStr = nameStr;
        end
    end

    methods(Static)
        function [nameStrs, m] = getNameStrs()
            [m,~] = enumeration('LaunchVehicleExtremaTypeEnum');

            nameStrs = {};
            for(i=1:length(m)) %#ok<*NO4LP>
                nameStrs{end+1} = m(i).nameStr; %#ok<AGROW>
            end
        end

        function [ind, mI] = getIndOfListboxStrs(recordingEnum)
            [m,~] = enumeration('LaunchVehicleExtremaTypeEnum');

            ind = -1;
            mI = [];
            for(i=1:length(m))
                if(m(i) == recordingEnum)
                    ind = i;
                    mI = m(i);
                    break;
                end
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef Pooling < dagnn.Filter
  properties
    method = 'max'
    poolSize = [1 2]
    opts = {'cuDNN'}
  end

  methods
    function outputs = forward(self, inputs, params)
      outputs{1} = vl_nnpool(inputs{1}, self.poolSize, ...
                             'pad', self.pad, ...
                             'stride', self.stride, ...
                             'method', self.method, ...
                             self.opts{:}) ;
    end

    function [derInputs, derParams] = backward(self, inputs, params, derOutputs)
      derInputs{1} = vl_nnpool(inputs{1}, self.poolSize, derOutputs{1}, ...
                               'pad', self.pad, ...
                               'stride', self.stride, ...
                               'method', self.method, ...
                               self.opts{:}) ;
      derParams = {} ;
    end

    function kernelSize = getKernelSize(obj)
      kernelSize = obj.poolSize ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes = getOutputSizes@dagnn.Filter(obj, inputSizes) ;
      outputSizes{1}(3) = inputSizes{1}(3) ;
    end

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

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function arccosh_values_test ( )

%*****************************************************************************80
%
%% ARCCOSH_VALUES_TEST demonstrates the use of ARCCOSH_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ARCCOSH_VALUES_TEST:\n' );
  fprintf ( 1, '  ARCCOSH_VALUES stores values of \n' );
  fprintf ( 1, '  the hyperbolic arc cosine function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = arccosh_values ( n_data );

    if ( n_data == 0 )
      break
    end

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

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%kkappend 'Append (or Prepend) Input Objects to Base Object'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kappend.pane file
%
% Parameters:
% InputFile: i1 'Input 1 (Base)', required: 'Base Input data object (input 1)'
% InputFile: i2 'Input 2', optional: 'input file 2'
% InputFile: i3 'Input 3', optional: 'input file 3'
% InputFile: i4 'Input 4', optional: 'input file 4'
% InputFile: i5 'Input 5', optional: 'input file 5'
% InputFile: i6 'Input 6', optional: 'input file 6'
% InputFile: i7 'Input 7', optional: 'input file 7'
% InputFile: i8 'Input 8', optional: 'input file 8'
% InputFile: i9 'Input 9', optional: 'input file 9'
% InputFile: i10 'Input 10', optional: 'input file 10'
% OutputFile: o 'Output', required: 'Resulting output data object'
% Toggle: w 'Width', default: 0: 'Stack objects along width dimension'
% Toggle: h 'Height', default: 0: 'Stack objects along height dimension'
% Toggle: d 'Depth', default: 0: 'Stack objects along depth dimension'
% Toggle: t 'Time', default: 0: 'Stack objects along time dimension'
% Toggle: e 'Elements', default: 0: 'Stack objects along elements dimension'
%
% Example: o = kkappend({i1, i2, i3, i4, i5, i6, i7, i8, i9, i10}, {'i1','';'i2','';'i3','';'i4','';'i5','';'i6','';'i7','';'i8','';'i9','';'i10','';'o','';'w',0;'h',0;'d',0;'t',0;'e',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kappend - Append (or Prepend) Input Objects to Base Object
%
%  DESCRIPTION
% The "Append" operator (kappend) combines all input data objects
% (i1 to i10) into a single output object (o).
% The objects are combined in order from Input 1 to Input 10.  If the
% user specifies placement as "Prepend Objects to Base Object" (-append 0),
% then the first input will be the last appended to the output object.
% The dimension flags Width (w_flag), Height (h_flag), Depth (d_flag),
% Time (t_flag), and Elements (e_flag), specify which dimension the objects
% will be stacked along.  For example, stacking two images along width
% would produce an output with the images side by side.  One and only one
% of the dimension flags can be true.
%
%  "Value Data" 5
% All input objects must have value data.  The value size of each input
% object will be zero-padded to the maximum of all sizes, except in the
% stacking dimension.  The output size in the stacking dimension will
% be the sum of that dimension across all input objects.  The output object
% data type will be the highest data type of the input data objects.
%
%  "Map Data" 5
% When some or all of the input objects have map data, the treatment of the
% maps, and how the data is represented in the output object, depends on the
% mapping option (mapping) specified by the user.
% Possible mapping options are:  (0) Map Data Thru Maps, (1) Use First Map Only,
% and (2) Create a Map for Each Input.
% If no input objects have maps, the mapping option is ignored.
% If there are doubts about which
% mapping option to use, the safest bet is to map the data thru the maps.
%
%
%  "Map Data Thru Maps:" 5
% 	All data will be mapped before the data objects are combined.  The output
% 	object will not have a map.
%
%  "Use First Map Only:" 5
% 	In this case, the map data and color attributes of the first input
% 	object that has a map are directly transferred to the output object.
% 	The depth, time, and elements dimensions of that map must be 1, otherwise
% 	this mapping option will fail.  Note that by selecting this mapping
% 	option, you are assuming that the value segments of all objects being
% 	appended have valid indices into that map.
%
%  "Create a Map for Each Input:" 5
% 	(All input objects may have maps.)
% 	For each input object that does not have a map, one will be created for it.
% 	The the output object color attributes will be copied from the first
% 	input object that has a map.  New map values will be assigned using the
% 	autocolor attribute.
%
%
% 	The map dimension corresponding to the stacking dimension (Depth, Time,
% 	or Elements) will reflect the new size of the value data in that dimension,
% 	and the maps from each object will be appended accordingly.
% 	The maximum width and height of all maps will be propagated to the output,
% 	with zero-padding being applied where necessary.
% 	Padding in the other map dimensions may occur as well.
% 	For example, if objects with the following map and value
% 	Depth sizes are being stacked along the Depth dimension,
%
% 			obj1  DEP=3  MAP_DEP = 3
% 			obj2  DEP=4  MAP_DEP = 1
% 			obj2  DEP=2  MAP_DEP = 2
%
% 	the resulting value data will have a depth of 9, so the map depth will be 9
% 	too.  In this example, the second object's map will be zero-padded to have
% 	a map depth of four.
%
%
% 	If the objects are being combined along the width or
% 	height dimensions, this mapping option will fail.
%
%
%
% The output object map data type will be the highest of the input object map
% data types.
%
%  "Mask Data" 5
% If any input object has mask data, then the output
% object will also have mask data.  The mask data will be resized
% in the same manner as value data is resized.  The mask will be padded
% as valid (1). The highest
% data type of input masks will be propagated to the output mask.
%
%  "Location and Time Data" 5
% In all cases, the time and location data and/or attributes
% of the first input object will transferred to the output
% object.  If location or time data is propagated, data values
% along the stacking dimension will not be changed, and
% therefore, the location or time data may become invalid.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::kinset
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkappend(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkappend(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'i3', '__input';'i4', '__input';'i5', '__input';'i6', '__input';'i7', '__input';'i8', '__input';'i9', '__input';'i10', '__input';'o', '__output';'w', 0;'h', 0;'d', 0;'t', 0;'e', 0};
maxval={0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0};
minval={0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0};
istoggle=[0,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','InputFile','InputFile','InputFile','InputFile','InputFile','InputFile','InputFile','InputFile','OutputFile','Toggle','Toggle','Toggle','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kappend"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
addpath('../toolbox/');
addpath('../image_blur/');
addpath('../blur_functions//');
addpath('../convolutional_wasserstein/');
addpath('../../data/images/shapes/'); % low-res images

rep = '../results/shapes/';
if not(exist(rep))
    mkdir(rep);
end

%% Read data

if not(exist('imresize'))
    imresize = @(x,s)image_resize(x,s);
end

targetSize = 199;

shape_list = [2 3 4];
nImages = length(shape_list);

data = zeros(targetSize*targetSize,nImages);
for i=1:nImages
    im = imread(sprintf('shape%dfilled.png', shape_list(i)));
    im = im2double(rgb2gray(im));

    padding = abs(size(im,2)-size(im,1));
    pad1 = floor(padding/2);
    pad2 = padding-pad1;

    if size(im,1) < size(im,2)
        im = [zeros(pad1,size(im,2)) ; im ; zeros(pad2,size(im,2))];
    elseif size(im,2) < size(im,1)
        im = [zeros(size(im,1),pad1) im zeros(size(im,1),pad2)];
    end

    im = 1-imresize(im,[targetSize targetSize]);
    im(im<0) = 0;

    im = im > .01;

    im = im + 1e-3;
    im = im / sum(im(:));

    % figure;imagesc(-im);axis equal;axis off;colormap gray;

    data(:,i) = im(:);
end

n = size(data,1);
areaWeights = ones(n,1)/n;
data = data*n;

entropies = -sum(bsxfun(@times,data.*log(data),areaWeights),1);
maxEntropy = max(entropies);

%% Set up blur

filterSize = 1.5; %was 5
imSize = [targetSize targetSize];

if exist('imfilter')
    % using image toolbox
    h = fspecial('gaussian',[1 max(imSize)],filterSize);% hsize sigma
    h = h / sum(h);
    imBlur = @(x) imfilter(imfilter(x,h,'replicate'),h','replicate');
else
    blur = load_filtering('imgaussian', targetSize);
    imBlur = @(x)blur(x,filterSize);
end

imagesc(imBlur(reshape(data(:,1),imSize)));

blurColumn = @(x) reshape(imBlur(reshape(x,imSize)),[],1);
blurAll2 = @(x) cell2mat(cellfun(blurColumn, num2cell(x,1), 'UniformOutput', false));
blurAll = @(x) blurAll2(blurAll2(x));

% blurColumn = @(x) reshape(fastBlur(reshape(x,[targetSize,targetSize]),filterSize),[],1);
% blurAll = @(x) cell2mat(cellfun(blurColumn, num2cell(x,1), 'UniformOutput', false));

imagesc(reshape(blurAll(data(:,1)),imSize));
axis equal;
axis off;

%% Compute barycenter

entropyLimit = Inf;

steps = linspace(0,1,5);% linspace(-.5,1.5,9);

averages = cell(length(steps),length(steps));
barycenters = cell(length(steps),length(steps));

options.niter = 400; %not 300
options.verb = 2;
options.tol = 1e-15;
resh = @(x)reshape(x, imSize);
options.disp = @(x)imageplot(-resh(x));

W = { ...
    [0, 0, 1] ...
    [1, 0, 3] [0, 1, 3] ...
    [1,0,1] [1,1,2] [0,1,1] ...
    [3,0,1] [2,1,1] [1,2,1] [0,3,1] ...
    [1,0,0] [3,1,0] [1,1,0] [1,3,0] [0,1,0] ...
    };

B = {};
for i=1:length(W)
    alpha = W{i}; alpha = alpha/sum(alpha);
    B{i} = convolutionalBarycenter(data,alpha,areaWeights,blurAll,blurAll,entropyLimit,options);
    u = rescale(-resh(B{i}));
    imwrite(u, [rep 'shape-barycenter-' num2str(i) '.png'], 'png');
    imwrite(double(u>.5), [rep 'tresh-barycenter-' num2str(i) '.png'], 'png');
end

author: gpeyre


text:
function disk_integrals_test01 ( )

%*****************************************************************************80
%
%% DISK_INTEGRALS_TEST01 uses DISK01_SAMPLE to estimate various integrals.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2014
%
%  Author:
%
%    John Burkardt
%
  m = 2;
  n = 4192;
  test_num = 20;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  Estimate monomial integrals using Monte Carlo\n' );
  fprintf ( 1, '  over the interior of the unit disk in 2D.\n' );
%
%  Get sample points.
%
  seed = 123456789;
  [ x, seed ] = disk01_sample ( n, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of sample points used is %d\n', n );
%
%  Randomly choose X,Y exponents between 0 and 8.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  If any exponent is odd, the integral is zero.\n' );
  fprintf ( 1, '  We will restrict this test to randomly chosen even exponents.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ex  Ey     MC-Estimate           Exact      Error\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    [ e, seed ] = i4vec_uniform_ab ( m, 0, 4, seed );

    e(1:m) = e(1:m) * 2;

    value = monomial_value ( m, n, e, x );

    result = disk01_area ( ) * sum ( value(1:n) ) / n;
    exact = disk01_monomial_integral ( e );
    error = abs ( result - exact );

    fprintf ( 1, '  %2d  %2d  %14.6g  %14.6g  %10.2e\n', ...
      e(1:m), result, exact, error );

  end

  return
end


text:
%  Figure 10.17      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
%  fig10_17.m is a script to generate Fig. 10.17, transient
%  response of the PD plus notch compensator of the
%  satellite position control, non-colocated case with the stiff-spring
%  parameters
m=[1 .1] ; k1 = [0 .4] ; b1=[0 .04*sqrt(.04)];
[f,g,h,j]=twomass(m,k1,b1);

no3=conv(.25*[2 1],[1/.81 0 1]);
do3=conv([1/40 1],[1/625 2/25 1]);
[Ac3,Bc3,Cc3,Dc3]=tf2ss(no3,do3);
sys1=ss(f,g,h,j);
sys3=ss(Ac3,Bc3,Cc3,Dc3);
sysol=series(sys3,sys1);
[Aol,Bol,Col,Dol]=ssdata(sysol);
[Acl]=[Aol-Bol*Col];
hold off; clf
%subplot(224);
t=0:0.1:40;
sys=ss(Acl,Bol,Col,Dol);
[y]=step(sys,t);
plot(t,y);
xlabel('Time (sec)');
ylabel('Amplitude');
title('Fig. 10.17 Closed-loop step response for KD_3(s)Ghat(s)')
grid;

author: Sable

repo: mcbench-benchmarks


text:
function obj = Make_f15( obj, ts, te, dt, varargin )
% obj = Make_f15(obj,ts, te, dt, varargin)
% Input a msh class object, the simulation start time, end times, timestep
% and any other optional arguments and get back a mesh class object with
% the f15 structure populated
%
%  Inputs:   1) A msh class obj (requires nodestrings if you want to get
%               the elevation specified boundary forcings), requires f13
%               populated to put in the NWP names if you want any f13 options.
%            2) ts: the simulation start time in datestring format e.g.
%               ts = '01-Jan-2018 03:00'
%            3) te: the simulation end time in datestring format e.g.
%               te = '01-Feb-2018 22:00'
%            4) dt: the simulation time step in seconds
%            5) Optional name-value arguments:
%
%               'const' : the harmonic constituents to use for the tidal
%               potential, elevation specified boundaries (if you have
%               nodestrings and a tidal database specified), and the
%               harmonic analysis. Give value as a cell array of characters
%               e.g: C = {'M2','S2','K1','O1'}; or a string array,
%               e.g: C = ["M2","S2","K1","O1"]; or u can put:
%                    C = "major8" to do all the major eight constituents
%               Note that the tidal nodal factors and so forth will be
%               automatically calculated based on ts and te.
%
%               'tidal_database' : the directory + filename tidal database of
%               elevations to interpolate to your boundaries.
%               This option requires the 'const' name-value pair input.
%               Use either:
%               a) TPXO9.1 (Netcdf), available for download at:
%                  http://volkov.oce.orst.edu/tides/global.html. The
%                  filename will be 'h_tpxo9.v1.nc'.
%               b) The TPXO Atlas (Netcdf) in which case the filename
%                  is 'h_**_tpxo9_atlas_30.nc' (for e.g., the tpxo9 atlas)
%                  where the wildcard is used in place of the constituent name
%
%               'sta_database' : a cell, where the first element contains
%               a character referring to a specific database. At the moment
%               it only handles CO-OPS/NOS/NOAA stations. The second
%               element contains a scalar either 1,2,3 corresponding to
%               elev, vel and met stations. e.g. set
%               ...('sta_database',{'CO-OPS',1,'NDBC',2,'NDBC',3}) to output
%               the CO-OPS stations for elevation and NDBC stations for velocity and
%               meteorological recording
%
%  Outputs:  1) msh class obj with f15 struct populated
%
%  Author:      William Pringle and Keith Roberts
%  Created:     March 15 2018
%  Updated:     July 15, 2021 by Shintaro Bunya and Keith Roberts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Putting in default options if f15 isn't already populated
if isempty(obj.f15)
    % Run Description and ID
    f15dat.rundes = obj.title ;  % RUNDES
    f15dat.runid  = 'Run_001' ;  % RUNID
    % NFOVER
    f15dat.nfover = 1 ;
    % NABOUT
    f15dat.nabout = 0 ;
    % NSCREEN
    f15dat.nscreen = 0;
    % IHOT (cold start)
    f15dat.ihot = 0;
    % ICS
    if all(abs(obj.p(:,1))) <= 180 && all(abs(obj.p(:,2))) <= 90
        f15dat.ics = 2 ; % uses lat, lon (probably)
    else
        f15dat.ics = 1 ; % not lat, lon
    end
    % IM (default is explicit with smagorinsky)
    f15dat.im = 511112 ;
    f15dat.iden = [];
    % NOLIBF (normal bottom friction)
    f15dat.nolibf = 1 ;
    % NOLIFA (fininte amplitude terms and wetting-drying on)
    f15dat.nolifa = 2 ;
    % NOLICA (the advection terms on)
    f15dat.nolica = 1 ;
    % NOLICAT (the advection terms on)
    f15dat.nolicat = 1 ;
    % NWP
    f15dat.nwp = 0 ;
    % NCOR
    if f15dat.ics == 2
        f15dat.ncor = 1 ; %spatially varying coriolis from latitude
    else
        f15dat.ncor = 0 ; % read in Coriolis
    end
    % NTIP
    f15dat.ntip = 0 ;
    % NWS
    f15dat.nws = 0 ;
    f15dat.wtimnc = 0;
    % NRAMP
    f15dat.nramp = 0 ;
    % G
    f15dat.gravity = 9.81 ;
    % TAU0 (the automatic spatially varying tau0)
    f15dat.tau0 = -1 ;
    f15dat.tau0minmax = 0;
    % DTDP
    f15dat.dtdp = 0 ;
    % STATIM
    f15dat.statim = 0 ;
    % REFTIM
    f15dat.reftim = 0 ;
    % RNDY
    f15dat.rndy = 0;
    % DRAMP
    f15dat.dramp = 0;
    % A00, B00, C00
    f15dat.a00b00c00 = [0 1 0] ;
    % H0
    f15dat.h0 = [0.1 0 0 0.01];
    % SLAM0, SFEA0
    f15dat.slam = mean(obj.p) ;
    % CF
    f15dat.taucf = 0.0025 ;
    % ESLM (smag 0.05 coefficient)
    f15dat.elsm = -0.05;
    % CORI
    f15dat.cori = 0;
    % NTIF
    f15dat.ntif = 0;
    % NBFR
    f15dat.nbfr = 0 ;
    %  ANGINN
    f15dat.anginn = 90 ;
    % Normal flow external boundary condition (needs to be fixed)
    ibtype = [2 12 22 32 52] ;
    sm = 0 ;
    if ~isempty(obj.bd)
        for k = 1: length(ibtype)
            sm = sm + sum(~(obj.bd.ibtype - ibtype(k))) ;
        end
    end
    if ( sm > 0 )
        f15dat.nffr = 0 ;
    end
    % NOUTE, TOUTSE, TOUTFE, NSPOOLE
    f15dat.oute = zeros(1,4) ;
    f15dat.elvstaloc = [];
    % NSTAE
    f15dat.nstae = 0 ;
    % NOUTV, TOUTV, TOUTFV, NSPOOLV
    f15dat.outv = zeros(1,4);
    f15dat.velstaloc = [];
    % NSTAV
    f15dat.nstav = 0 ;
    % NSTAM, TOUSM, TOUFM, NSPOOLM
    f15dat.outm = zeros(1,4) ;
    f15dat.metstaloc = [];
    % NSTAM
    f15dat.nstam = 0 ;

    % NOUTGE
    f15dat.outge = zeros(1,4) ;
    % NOUTGV
    f15dat.outgv = zeros(1,4) ;
    % NOUTGM
    f15dat.outgm = zeros(1,4);
    % NFREQ
    f15dat.nfreq = 0 ;
    % THAS, THAF, NHAINC, FMV
    f15dat.outhar = zeros(1,4) ;
    % NHASE, NHASV, NHAGE, NHAGV
    f15dat.outhar_flag = zeros(1,4) ;
    % NHSTAR, NHSINC
    f15dat.nhstar = zeros(1,2) ;
    % ITITER, ISLDIA, CONVCR, ITMAX
    f15dat.ititer = [1 0 1e-10 25] ;
    % Extra lines for NETCDF
    f15dat.nextraline = 10;
    f15dat.extraline(1).msg = obj.title ;
    f15dat.extraline(2).msg = 'Affiliation';
    f15dat.extraline(3).msg = 'Made by OceanMesh2D';
    f15dat.extraline(4).msg = 'History: None';
    f15dat.extraline(5).msg = 'https://github.com/CHLNDDEV/OceanMesh2D/';
    f15dat.extraline(6).msg = 'Comments: None' ;
    f15dat.extraline(7).msg = 'Host: Name';
    f15dat.extraline(8).msg = 'Metric, Datum';
    f15dat.extraline(9).msg = 'name@instit.edu';
    f15dat.extraline(10).msg = '';

    % control lists
    % met
    f15dat.controllist(1).type = 'met';
    f15dat.controllist(1).var(1).name = 'WindDragLimit';
    f15dat.controllist(1).var(1).val = 0.0025;
    f15dat.controllist(1).var(2).name = 'PRBCKGRND';
    f15dat.controllist(1).var(2).val = 1013;
    f15dat.controllist(1).var(3).name = 'DragLawString';
    f15dat.controllist(1).var(3).val = 'default';
    f15dat.controllist(1).var(4).name = 'outputWindDrag';
    f15dat.controllist(1).var(4).val = 'F';
    f15dat.controllist(1).var(5).name = 'invertedBarometerOnElevationBoundary';
    f15dat.controllist(1).var(5).val = 'F';
    % dynamicwaterlevelcorrection control
    f15dat.controllist(2).type = 'dynamicWaterLevelCorrection';
    f15dat.controllist(2).var(1).name = [f15dat.controllist(2).type 'FileName'];
    f15dat.controllist(2).var(1).val = 'null';
    f15dat.controllist(2).var(2).name = [f15dat.controllist(2).type 'Multiplier'];
    f15dat.controllist(2).var(2).val = 1.0;
    f15dat.controllist(2).var(3).name = [f15dat.controllist(2).type 'RampStart'];
    f15dat.controllist(2).var(3).val = 0.0;
    f15dat.controllist(2).var(4).name = [f15dat.controllist(2).type 'RampEnd'];
    f15dat.controllist(2).var(4).val = 0.0;
    f15dat.controllist(2).var(5).name = [f15dat.controllist(2).type 'RampReferenceTime'];
    f15dat.controllist(2).var(5).val = 'coldstart';
    f15dat.controllist(2).var(6).name = [f15dat.controllist(2).type 'SkipSnaps'];
    f15dat.controllist(2).var(6).val = 0;

    % Put into the msh class
    obj.f15 = f15dat;
end

%% Changing following parameters based on standard require inputs and
%% populated properties of the msh class
% NSCREEN (once per simulation hour)
obj.f15.nscreen = round(3600/dt);
% DTDP (time step)
obj.f15.dtdp = dt ;
% Runday
obj.f15.rndy = datenum(te)-datenum(ts);

obj.f15.extraline(10).msg = datestr(ts, 'yyyy-mm-dd HH:MM:SS UTC');

% Checking f11
if ~isempty(obj.f11)
    % using WJP's new IM for a 2DDI Baroclinic run
    obj.f15.im = 511113 ;
    obj.f15.iden = 1;
end

% Checking f13
if ~isempty(obj.f13)
    obj.f15.nwp = obj.f13.nAttr ;
    if ( obj.f15.nwp > 0 )
        for l = 1: obj.f15.nwp
            obj.f15.AttrName(l).name = obj.f13.defval.Atr(l).AttrName ;
            if strfind(obj.f15.AttrName(l).name,'primitive')
                obj.f15.tau0 = -3; % change tau0 option
            end
        end
    end
end

%% Test optional arguments
% default
const = [];
tidal_database = [];
sta_database = [];
if ~isempty(varargin)
    names = {'const','tidal_database','NWS','WTIMNC','tau0','tau0minmax',...
             'sta_database'};
    for ii = 1:length(names)
        ind = find(~cellfun(@isempty,strfind(varargin(1:2:end),names{ii})));
        if ~isempty(ind)
            if ii == 1
                const = varargin{ind*2};
            elseif ii == 2
                tidal_database = varargin{ind*2};
            elseif ii == 3
                obj.f15.nws = varargin{ind*2};
            elseif ii == 4
                obj.f15.wtimnc = varargin{ind*2};
            elseif ii == 5
                obj.f15.tau0 = varargin{ind*2};
            elseif ii == 6
                obj.f15.tau0minmax = varargin{ind*2};
            elseif ii == 7
                sta_database = varargin{ind*2};
            end
        end
    end
end

%% Changing following parameters based on optional argument inputs
% WTIMINC
if obj.f15.nws > 0 %~mod(obj.f15.nws - 1,10) ||  ~mod(obj.f15.nws - 11,10) || obj.f15.nws == 8
    if obj.f15.wtimnc == 0
        error(['Must specify WTIMNC option with NWS = ' num2str(obj.f15.nws)])
    end
end

% Tau0FullDomainMin, Tau0FullDomainMax
if ( abs(obj.f15.tau0 + 5.D0) < 1e-10 )
    if obj.f15.tau0minmax == 0
        error(['Must specify tau0minmax with tau0 option = ' num2str(obj.f15.tau0)])
    end
end

% Tidal potential and harmonic analysis stuff
if ~isempty(const)
    % Get the tidal factors etc based on ts and te
    if ischar(const)
        const = {const};
    end
    obj = tide_fac(obj,ts,te,const);

    % Harmonic analysis stuff (copy in tidal potential)
    obj.f15.nfreq = obj.f15.ntif; K = [];
    for k = 1: obj.f15.nfreq
        obj.f15.harfreq(k).name = obj.f15.tipotag(k).name;
        % get only the frequencies, nodal factor and equilibrium argument
        obj.f15.harfreq(k).val = obj.f15.tipotag(k).val([2,4,5]);
        if isnan(obj.f15.tipotag(k).val(1))
            % get rid of this constituent from potential
            obj.f15.ntif = obj.f15.ntif - 1;
            K(end+1) = k;
        end
    end
    obj.f15.tipotag(K) = [];

    if obj.f15.ntip == 0
        disp('Setting ntip = 1')
        obj.f15.ntip = 1;
    end
end

% Elevation Specified Boundary Conditions (tides)
if ~isempty(const) && ~isempty(tidal_database)
    % copy in the harmonic analysis stuff
    obj.f15.nbfr = obj.f15.nfreq;
    for k = 1: obj.f15.nbfr
        obj.f15.bountag(k).name = obj.f15.harfreq(k).name;
        % get only the frequencies, nodal factor and equilibrium argument
        obj.f15.bountag(k).val = obj.f15.harfreq(k).val;
    end

    % Do the interpolation to the boundaries
    obj = tidal_data_to_ob(obj,tidal_database,{obj.f15.bountag.name});
end

if ~isempty(const)
   % Here for the SA tidal potential the equilibrium argument is now
   % corrected to be referenced to the vernal equinox (subtract of about
   % 77 deg). But the equilibrim argument for the boundary and the analysis
   % is kept as referenced to the Calender year.
   ii = find(startsWith({obj.f15.tipotag.name},'SA'));
   if ~isempty(ii)
       % need to refer to vernal equinox instead of calender year,
       % at least for the tidal potential (subtract about 76.9 deg);
       obj.f15.tipotag(ii).val(5) = obj.f15.tipotag(ii).val(5) - 78/365*360;
       if obj.f15.tipotag(ii).val(5) < 0
          obj.f15.tipotag(ii).val(5) = obj.f15.tipotag(ii).val(5) + 360;
       end
   end
end
%
if ~isempty(sta_database)
    % parse the station reader
    obj = tidal_stations_parser(obj,string(sta_database(1:2:end-1)),...
                                sta_database(2:2:end));
end
%
% % NOUTGC
% if ( f15dat.im == 10 )
%     f15dat.outgc = readlinevec( fid ) ;
% end
%
% % NOUTGW
% if ( f15dat.nws ~= 0 )
%     f15dat.outgw = readlinevec( fid ) ;
% end

end

author: CHLNDDEV


text:
function chebyshev_polynomial_test15 ( )

%*****************************************************************************80
%
%% CHEBYSHEV_POLYNOMIAL_TEST15 tests V_POLYNOMIAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBYSHEV_POLYNOMIAL_TEST15:\n' );
  fprintf ( 1, '  V_POLYNOMIAL_VALUES stores values of\n' );
  fprintf ( 1, '  the Chebyshev polynomials.\n' );
  fprintf ( 1, '  V_POLYNOMIAL evaluates the polynomial.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                        Tabulated                 Computed\n' );
  fprintf ( 1, '     N        X           V(n,x)                    V(n,x)                     Error\n' );

  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx1 ] = v_polynomial_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2_vec = v_polynomial ( 1, n, x );
    fx2 = fx2_vec(1,n+1);
    e = fx1 - fx2;

    fprintf ( 1, '  %4d  %12f  %24.16e  %24.16e  %8.2g\n', n, x, fx1, fx2, e );

  end

  return
end


text:
%-----------------------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 4899 $)
%-----------------------------------------------------------------------
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_entry.type = 'cfg_entry';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_entry.name = 'String';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_entry.tag = 'instr';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_entry.strtype = 's';
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_entry.extras = [];
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_entry.num = [0 Inf];
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_entry.check = [];
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_entry.help = {'Enter a string.'};
matlabbatch{1}.menu_cfg{1}.menu_entry{1}.conf_entry.def = [];
matlabbatch{2}.menu_cfg{1}.menu_entry{1}.conf_entry.type = 'cfg_entry';
matlabbatch{2}.menu_cfg{1}.menu_entry{1}.conf_entry.name = 'Evaluated Input';
matlabbatch{2}.menu_cfg{1}.menu_entry{1}.conf_entry.tag = 'ineval';
matlabbatch{2}.menu_cfg{1}.menu_entry{1}.conf_entry.strtype = 'e';
matlabbatch{2}.menu_cfg{1}.menu_entry{1}.conf_entry.extras = [];
matlabbatch{2}.menu_cfg{1}.menu_entry{1}.conf_entry.num = [];
matlabbatch{2}.menu_cfg{1}.menu_entry{1}.conf_entry.check = [];
matlabbatch{2}.menu_cfg{1}.menu_entry{1}.conf_entry.help = {'Enter an evaluated input.'};
matlabbatch{2}.menu_cfg{1}.menu_entry{1}.conf_entry.def = [];
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.type = 'cfg_files';
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.name = 'NIfTI Images';
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.tag = 'innifti';
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.filter = 'image';
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.ufilter = '.*';
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.dir = '';
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.num = [0 Inf];
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.check = [];
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.help = {'Select NIfTI Images'};
matlabbatch{3}.menu_cfg{1}.menu_entry{1}.conf_files.def = [];
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.type = 'cfg_files';
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.name = 'MATLAB .mat Files';
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.tag = 'inmat';
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.filter = 'mat';
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.ufilter = '.*';
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.dir = '';
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.num = [0 Inf];
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.check = [];
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.help = {'Select MATLAB .mat files.'};
matlabbatch{4}.menu_cfg{1}.menu_entry{1}.conf_files.def = [];
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.type = 'cfg_files';
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.name = 'Any Files';
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.tag = 'inany';
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.filter = 'any';
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.ufilter = '.*';
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.dir = '';
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.num = [0 Inf];
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.check = [];
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.help = {'Select any kind of files.'};
matlabbatch{5}.menu_cfg{1}.menu_entry{1}.conf_files.def = [];
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.type = 'cfg_files';
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.name = 'Directory';
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.tag = 'indir';
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.filter = 'dir';
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.ufilter = '.*';
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.dir = '';
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.num = [1 Inf];
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.check = [];
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.help = {'Directory'};
matlabbatch{6}.menu_cfg{1}.menu_entry{1}.conf_files.def = [];
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.type = 'cfg_repeat';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.name = 'Job Inputs';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.tag = 'inputs';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1) = cfg_dep;
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).tname = 'Values Item';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).tgt_spec = {};
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).sname = 'Entry: String (cfg_entry)';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).src_exbranch = substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).src_output = substruct('()',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{2}(1) = cfg_dep;
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{2}(1).tname = 'Values Item';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{2}(1).tgt_spec = {};
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{2}(1).sname = 'Entry: Evaluated Input (cfg_entry)';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{2}(1).src_exbranch = substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{2}(1).src_output = substruct('()',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{3}(1) = cfg_dep;
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{3}(1).tname = 'Values Item';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{3}(1).tgt_spec = {};
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{3}(1).sname = 'Files: NIfTI Images (cfg_files)';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{3}(1).src_exbranch = substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{3}(1).src_output = substruct('()',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{4}(1) = cfg_dep;
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{4}(1).tname = 'Values Item';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{4}(1).tgt_spec = {};
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{4}(1).sname = 'Files: MATLAB .mat Files (cfg_files)';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{4}(1).src_exbranch = substruct('.','val', '{}',{4}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{4}(1).src_output = substruct('()',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{5}(1) = cfg_dep;
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{5}(1).tname = 'Values Item';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{5}(1).tgt_spec = {};
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{5}(1).sname = 'Files: Any Files (cfg_files)';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{5}(1).src_exbranch = substruct('.','val', '{}',{5}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{5}(1).src_output = substruct('()',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{6}(1) = cfg_dep;
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{6}(1).tname = 'Values Item';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{6}(1).tgt_spec = {};
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{6}(1).sname = 'Files: Directory (cfg_files)';
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{6}(1).src_exbranch = substruct('.','val', '{}',{6}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{6}(1).src_output = substruct('()',{1});
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.num = [0 Inf];
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.forcestruct = false;
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.check = [];
matlabbatch{7}.menu_cfg{1}.menu_struct{1}.conf_repeat.help = {'Assemble the set of input items for one run of the job.'};
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.type = 'cfg_repeat';
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.name = 'Runs';
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.tag = 'runs';
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1) = cfg_dep;
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).tname = 'Values Item';
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).tgt_spec = {};
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).sname = 'Repeat: Job Inputs (cfg_repeat)';
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).src_exbranch = substruct('.','val', '{}',{7}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.values{1}(1).src_output = substruct('()',{1});
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.num = [1 Inf];
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.forcestruct = false;
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.check = [];
matlabbatch{8}.menu_cfg{1}.menu_struct{1}.conf_repeat.help = {'Repeat "Job Inputs" for each run of the job, even if you want to specify no inputs in this batch itself. The number of "Job Inputs" items determines how often this job is run.'};
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.type = 'cfg_files';
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.name = 'Job File(s)';
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.tag = 'jobs';
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.filter = 'batch';
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.ufilter = '.*';
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.dir = '';
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.num = [1 Inf];
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.check = [];
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.help = {'Select the job template(s). If multiple files are selected, they will be concatenated in selection order to form one job.'};
matlabbatch{9}.menu_cfg{1}.menu_entry{1}.conf_files.def = [];
matlabbatch{10}.menu_cfg{1}.menu_entry{1}.conf_const.type = 'cfg_const';
matlabbatch{10}.menu_cfg{1}.menu_entry{1}.conf_const.name = 'Don''t Save';
matlabbatch{10}.menu_cfg{1}.menu_entry{1}.conf_const.tag = 'dontsave';
matlabbatch{10}.menu_cfg{1}.menu_entry{1}.conf_const.val{1} = false;
matlabbatch{10}.menu_cfg{1}.menu_entry{1}.conf_const.check = [];
matlabbatch{10}.menu_cfg{1}.menu_entry{1}.conf_const.help = {'Do not save the generated jobs.'};
matlabbatch{10}.menu_cfg{1}.menu_entry{1}.conf_const.def = [];
matlabbatch{11}.menu_cfg{1}.menu_entry{1}.conf_entry.type = 'cfg_entry';
matlabbatch{11}.menu_cfg{1}.menu_entry{1}.conf_entry.name = 'Batch Filename Stub';
matlabbatch{11}.menu_cfg{1}.menu_entry{1}.conf_entry.tag = 'outstub';
matlabbatch{11}.menu_cfg{1}.menu_entry{1}.conf_entry.strtype = 's';
matlabbatch{11}.menu_cfg{1}.menu_entry{1}.conf_entry.extras = [];
matlabbatch{11}.menu_cfg{1}.menu_entry{1}.conf_entry.num = [1 Inf];
matlabbatch{11}.menu_cfg{1}.menu_entry{1}.conf_entry.check = [];
matlabbatch{11}.menu_cfg{1}.menu_entry{1}.conf_entry.help = {'The output filename will be generated by appending a job counter to this string.'};
matlabbatch{11}.menu_cfg{1}.menu_entry{1}.conf_entry.def = [];
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.type = 'cfg_files';
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.name = 'Batch Directory';
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.tag = 'outdir';
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.filter = 'dir';
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.ufilter = '.*';
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.dir = '';
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.num = [1 1];
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.check = [];
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.help = {'The generated batches will be saved into this folder.'};
matlabbatch{12}.menu_cfg{1}.menu_entry{1}.conf_files.def = [];
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.type = 'cfg_branch';
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.name = 'Save';
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.tag = 'savejobs';
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{1}(1) = cfg_dep;
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{1}(1).tname = 'Val Item';
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{1}(1).tgt_spec = {};
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{1}(1).sname = 'Batch Filename Stub (cfg_entry)';
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{1}(1).src_exbranch = substruct('.','val', '{}',{11}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{1}(1).src_output = substruct('()',{1});
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{2}(1) = cfg_dep;
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{2}(1).tname = 'Val Item';
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{2}(1).tgt_spec = {};
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{2}(1).sname = 'Batch Directory (cfg_files)';
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{2}(1).src_exbranch = substruct('.','val', '{}',{12}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.val{2}(1).src_output = substruct('()',{1});
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.check = [];
matlabbatch{13}.menu_cfg{1}.menu_struct{1}.conf_branch.help = {'Specify filename stub and output directory to save the generated files.'};
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.type = 'cfg_choice';
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.name = 'Save Generated Batch Jobs';
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.tag = 'save';
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{1}(1) = cfg_dep;
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{1}(1).tname = 'Values Item';
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{1}(1).tgt_spec = {};
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{1}(1).sname = 'Save (cfg_branch)';
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{1}(1).src_exbranch = substruct('.','val', '{}',{13}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{1}(1).src_output = substruct('()',{1});
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{2}(1) = cfg_dep;
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{2}(1).tname = 'Values Item';
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{2}(1).tgt_spec = {};
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{2}(1).sname = 'Don''t Save (cfg_const)';
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{2}(1).src_exbranch = substruct('.','val', '{}',{10}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.values{2}(1).src_output = substruct('()',{1});
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.check = [];
matlabbatch{14}.menu_cfg{1}.menu_struct{1}.conf_choice.help = {'The generated batch jobs can be saved for reference or debugging purposes.'};
matlabbatch{15}.menu_cfg{1}.menu_entry{1}.conf_menu.type = 'cfg_menu';
matlabbatch{15}.menu_cfg{1}.menu_entry{1}.conf_menu.name = 'Missing Inputs';
matlabbatch{15}.menu_cfg{1}.menu_entry{1}.conf_menu.tag = 'missing';
matlabbatch{15}.menu_cfg{1}.menu_entry{1}.conf_menu.labels = {
                                                              'Skip jobs with missing inputs, run filled jobs'
                                                              'Don''t run any jobs if missing inputs'
}';
matlabbatch{15}.menu_cfg{1}.menu_entry{1}.conf_menu.values = {
                                                              'skip'
                                                              'error'
}';
matlabbatch{15}.menu_cfg{1}.menu_entry{1}.conf_menu.check = [];
matlabbatch{15}.menu_cfg{1}.menu_entry{1}.conf_menu.help = {'Jobs with missing inputs (e.g. because of wrong input contents) can be skipped, while filled jobs are run. Alternatively, if any job has missing inputs, no jobs are run.'};
matlabbatch{15}.menu_cfg{1}.menu_entry{1}.conf_menu.def = [];
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.type = 'cfg_exbranch';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.name = 'Run Batch Jobs';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.tag = 'runjobs';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1) = cfg_dep;
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).tname = 'Val Item';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).tgt_spec = {};
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).sname = 'Job File(s) (cfg_files)';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).src_exbranch = substruct('.','val', '{}',{9}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{1}(1).src_output = substruct('()',{1});
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{2}(1) = cfg_dep;
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{2}(1).tname = 'Val Item';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{2}(1).tgt_spec = {};
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{2}(1).sname = 'Runs (cfg_repeat)';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{2}(1).src_exbranch = substruct('.','val', '{}',{8}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{2}(1).src_output = substruct('()',{1});
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{3}(1) = cfg_dep;
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{3}(1).tname = 'Val Item';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{3}(1).tgt_spec = {};
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{3}(1).sname = 'Save Generated Batch Jobs (cfg_choice)';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{3}(1).src_exbranch = substruct('.','val', '{}',{14}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{3}(1).src_output = substruct('()',{1});
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{4}(1) = cfg_dep;
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{4}(1).tname = 'Val Item';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{4}(1).tgt_spec = {};
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{4}(1).sname = 'Missing Inputs (cfg_menu)';
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{4}(1).src_exbranch = substruct('.','val', '{}',{15}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.val{4}(1).src_output = substruct('()',{1});
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.prog = @cfg_run_runjobs;
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.vout = @cfg_vout_runjobs;
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.check = [];
matlabbatch{16}.menu_cfg{1}.menu_struct{1}.conf_exbranch.help = {'Load a set of job files, fill missing inputs and run the filled job. This automates the creation and execution of batch jobs for a large number of identical computations.'};

author: spm


text:
function w = load_sac(request)
   %LOAD_SAC loads a waveform from SAC files
   % combineWaves isn't currently used!

   % request.combineWaves is ignored

   if isstruct(request)
      [thisSource, chanInfo, startTimes, endTimes, ~] = unpackDataRequest(request);
      for i=1:numel(chanInfo)
         for j=1:numel(startTimes)
            thisfilename = getfilename(thisSource,chanInfo(i),startTimes(j));
            w(i,j) = loadSacByFilenames(thisfilename);
            w(i,j) = extract(w(i,j), 'time', startTimes(j), endTimes(j)); % added by Glenn 2016/05/25
         end
      end
      w = w(:);

   else
      %request should be a filename
      % modified by Glenn 2016/05/25
      thisFilename = request;
      if exist(thisFilename, 'file')
          w = loadSacByFilenames(thisFilename);
      else
          w = waveform();
          warning(sprintf('File %s does not exist',thisFilename));
      end
   end

end

function w = loadSacByFilenames(filename)
   %LOADSAC  creates a waveform from a SAC file
   % waveform = loadSacByFilenames(waveform, filename);
   % To call, the first argument must be a waveform.  However this waveform
   % will never be accessed-- it will be overwritten by the SAC information.
   % Therefore, you can use "waveform" as your first argument, so that a
   % generic waveform object will be created for you
   %
   % If the filename is a cell of several SAC file names, such as:
   % fn = {'BOB.SHZ.AV','TOM.EHZ.AV','HARRY.EHZ.AV'}
   % then, w = loadSacByFilenames(waveform,fn) will return a matrix of waveforms, the
   % same size and shape  as fn.   fn(N,M) will be the SAC file name for the
   % waveform denoted by w(N,M).  In the example above, loadSacByFilenames will return a
   % 1x3 waveform object, with all three files.
   %
   % Programming suggestion???select files using a dialog box, by using the
   % following code:
   %   % display the dialog box used for selecting files, and
   %   % allow for multiple selections
   %   [f,d] = uigetfile({'*.*','all Files (*.*)'}, 'MultiSelect', 'on');
   %
   %   % f will have one or more file names, so concatenate them
   %   % with the directory
   %   filename = strcat(d,f);
   %
   %   %call loadSacByFilenames to load all these files into waveform objects
   %   w = loadSacByFilenames(waveform, filename);
   %
   % The SAC header is parsed out, with the following equivelent fields used
   % by waveform:
   % _waveform_ 	_sac_
   % STATION		KSTNM
   % CHANNEL		KCMPNM
   % FREQUENCY		1 / DELTA (ODELTA is not used)
   % START			NZYEAR, NZJDAY, NZHOUR, NZMIN, NZSEC, NZMSEC
   % UNITS			IDEP is parsed.  IUNK becomes "Counts"
   % All header fields with values - that is, header fields that do not
   % contain the value  -12345 (either numeric or text) - are put into
   % user-defined fields within the waveform.
   % These fields can then be accessed via get/set.  User-defined fields
   % includes those that are parsed into waveform's predefined  fields as
   % listed above.
   %
   % loadSacByFilenames uses routines originating in code borrowed from
   % Mike Thorn
   %
   % For more information about the SAC header, and the interpretation of the
   % fields, check out Lawrence Livermore's web site.  Search for "seismic
   % sac" and you'll find it.

   if ~iscell(filename)
      filename = {filename};
   end

   sacdata = readsacfile(filename{:});
   w = sac2waveform(sacdata);
   w = removeSacDuplicates(w);
   for n=1:numel(filename)
      w(n) = addhistory(w(n),'Loaded SAC file: %s',filename{n});
   end
end

function w = removeSacDuplicates(w)
   % because all the files are loaded with fields blindly filled in, several
   % of the duplicated (and easily outdated and contradictory) fields are
   % henceforth removed
   % These fields are regenerated when the SAC file is written back to disk.
   noHistory = true;
   w = delfield(w,'B',noHistory); %covered by get(w,'start_epoch');
   w = delfield(w,'E',noHistory); %covered by get(w,'end_epoch');
   w = delfield(w,'DEPMIN',noHistory); %covered by min(w);
   w = delfield(w,'DEPMAX',noHistory);%covered by max(w);
   w = delfield(w,'DEPMEN',noHistory);%covered by mean(w);
   w = delfield(w,'NPTS',noHistory); %covered by get(w,'data_length');
   w = delfield(w,'KSTNM',noHistory); %covered by get(w,'station');
   w = delfield(w,'KCMPNM',noHistory); %covered by get(w,'channel');
   w = delfield(w,'KNETWK',noHistory); %covered by get(w,'network');
   w = delfield(w,'DELTA',noHistory); %covered by get(w,'period');
   w = delfield(w,'NVHDR',noHistory); %covered by get(w,'data_length');
   w = delfield(w,'IDEP',noHistory); %covered by get(w,'units'); %idEP is just a code.

   w = delfield(w,'LEVEN',noHistory); %by nature, data is evenly spaced
   w = delfield(w,'IDEP',noHistory); %covered by get(w,'units'); %idEP is just a code.

   %the following date values are convenient, but also may be outdated easily.
   %w = delfield(w,'NZYEAR',noHistory); %covered by get(w,'start');
   %w = delfield(w,'NZJDAY',noHistory); %covered by get(w,'start');
   %w = delfield(w,'NZHOUR',noHistory); %covered by get(w,'start');
   %w = delfield(w,'NZMIN',noHistory); %covered by get(w,'start');
   %w = delfield(w,'NZSEC',noHistory); %covered by get(w,'start');
   %w = delfield(w,'NZMSEC',noHistory); %covered by get(w,'start');
end

author: geoscience-community-codes

repo: GISMO


text:
function [qrs,jpoints] = wqrsm_fast(data,fs,PWfreq,TmDEF,jflag)
% The swqrsm detector rewrited from swqrsm.c
% rewrited by Giulia Da Poian, August 7, 2018
% use struct swqrsm instead of global varaibles
% giulia.da.poian@emory.com
%
% input:
%   data:   ECG data, the swqrsm.c used the data in raw adus, if the input
%           data is in physical units, when the function found the peak-peak
%           value < 10, it will multiply the value by WFDB_DEFGAIN (200)
%   fs:     sampling frequency (default: 125)
%   PWfreq: power line (mains) frequency, in Hz (default: 60)
%   TmDEF:  minimum threshold value (default: 100)
%   jflag:  annotate J-points (ends of QRS complexes) (default: 0)
% output:
%   qrs:    QRS fiducial mark in samples
%   jpoints:J-points annotation, if jflag==1
%
% Matlab code based on:
%  file: swqrsm.c		Wei Zong      23 October 1998
% 			Last revised:  9 April 2010 (by G. Moody)
% -----------------------------------------------------------------------------
% swqrsm: Single-lead QRS detector based on length transform
% Copyright (C) 1998-2010 Wei Zong
%
% 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 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
% PARTICULAR PURPOSE.  See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place - Suite 330, Boston, MA 02111-1307, USA.
%
% You may contact the author by e-mail (wzong@mit.edu) or postal mail
% (MIT Room E25-505, Cambridge, MA 02139, USA).  For updates to this software,
% please visit PhysioNet (http://www.physionet.org/).
% ------------------------------------------------------------------------------
%
% This program analyzes an ECG signal, detecting QRS onsets and J-points, using
% a nonlinearly-scaled ECG curve length feature.  This version has been optimized
% for ECGs sampled at 125 Hz, but it can analyze ECGs sampled at any frequency
% using on-the-fly resampling provided by the WFDB library.
%
% `swqrsm' can process records containing any number of signals, but it uses only
% one signal for QRS detection (signal 0 by default;  this can be changed using
% the `-s' option, see below).  'swqrsm' has been optimized for adult human ECGs.
% For other ECGs, it may be necessary to experiment with the input sampling
% frequency and the time constants.
% ----------------------------------------------------------------------
% This matlab version has updates to deal with a variable sampling
% frequency but has not been exhanstively tested.
% ----- Add notes here:
%
%

if nargin<5
jflag = 0;
end
if nargin<4
TmDEF = 100;
end
if nargin<3
PWfreq = 60;
end
if nargin<2
fs = 125;
end

swqrsm.BUFLN = 16384; %	/* must be a power of 2, see ltsamp() */
EYE_CLS = 0.25; %   /* eye-closing period is set to 0.25 sec (250 ms) */
MaxQRSw = 0.13; %    /* maximum QRS width (130ms) */
NDP	= 2.5; %    /* adjust threshold if no QRS found in NDP seconds */
PWFreqDEF = PWfreq; %    /* power line (mains) frequency, in Hz (default) */
% TmDEF = 100; %	/* minimum threshold value (default) */
WFDB_DEFGAIN = 200.0; %  /* default value for gain (adu/physical unit) */

timer_d=0;
gain=WFDB_DEFGAIN;
swqrsm.lfsc = fix(1.25*gain*gain/fs);	% /* length function scale constant */
PWFreq = PWFreqDEF;	% /* power line (mains) frequency, in Hz */

% test data is physical units (mV) or raw adus units
datatest=data(1:fix(length(data)/fs)*fs);
if length(datatest)>fs
datatest=reshape(datatest,fs,[]);
end
test_ap=median(max(datatest)-min(datatest));
if test_ap<10 % peak-peak < 10 mV, may physical units
data=data*gain;
end

swqrsm.data = data;
swqrsm.lbuf = zeros(swqrsm.BUFLN,1);
swqrsm.ebuf = zeros(swqrsm.BUFLN,1);
swqrsm.ebuf(1:end)=fix(sqrt(swqrsm.lfsc));
swqrsm.lt_tt = 0;
swqrsm.aet = 0;
swqrsm.Yn = 0;
swqrsm.Yn1 = 0;
swqrsm.Yn2 = 0;

qrs = [];
jpoints = [];

Tm = fix(TmDEF / 5.0);
%     spm = 60 * fs;
%     next_minute = from + spm;
swqrsm.LPn = fix(fs/PWFreq); %		/* The LP filter will have a notch at the power line (mains) frequency */
if (swqrsm.LPn > 8)
swqrsm.LPn = 8;	% /* avoid filtering too agressively */
end
swqrsm.LP2n = 2 * swqrsm.LPn;
EyeClosing = fix(fs * EYE_CLS); % /* set eye-closing period */
ExpectPeriod = fix(fs * NDP);   % /* maximum expected RR interval */
swqrsm.LTwindow = fix(fs * MaxQRSw);   % /* length transform window size */

%     for i=1:2000
%         ltdata(i)=ltsamp(i);
%     end
%     /* Average the first 8 seconds of the length-transformed samples
%        to determine the initial thresholds Ta and T0. The number of samples
%        in the average is limited to half of the ltsamp buffer if the sampling
%        frequency exceeds about 2 KHz. */

t1 = fs*8;
if t1> fix(swqrsm.BUFLN*0.9)
    t1=swqrsm.BUFLN/2;
end

T0=0;
for t=1:t1
    swqrsm = ltsamp(t,swqrsm);
    T0 = T0 + swqrsm.lt_data;
end

T0 = T0 / t1;
Ta = 3 * T0;

% /* Main loop */
t=1;
learning=1;
while t<length(data)
    if (learning==1)
        if (t > t1)
            learning = 0;
            T1 = T0;
            t = 1;	% /* start over */
        else
            T1 = 2*T0;
        end
    end

    % /* Compare a length-transformed sample against T1. */
    swqrsm = ltsamp(t,swqrsm);
    if ( swqrsm.lt_data > T1) %	/* found a possible QRS near t */
        timer_d = 0; % /* used for counting the time after previous QRS */
        maxd = swqrsm.lt_data;
        mind = maxd;
        for tt = t+1:t + fix(EyeClosing/2)
            swqrsm = ltsamp(tt,swqrsm);
            if (swqrsm.lt_data > maxd)
                maxd = swqrsm.lt_data;
            end
        end
        for tt = t-1:-1:t - fix(EyeClosing/2)
            swqrsm = ltsamp(tt,swqrsm);
            if (swqrsm.lt_data < mind)
                mind = swqrsm.lt_data;
            end
        end
        if (maxd > mind+10) % /* There is a QRS near tt */
        % /* Find the QRS onset (PQ junction) */
            onset = fix(maxd/100) + 2;
            tpq = t - 5;
            for tt = t:-1:t - fix(EyeClosing/2)
                swqrsm = ltsamp(tt,swqrsm);
                swqrsm_1 = ltsamp(tt-1,swqrsm);
                swqrsm_2 = ltsamp(tt-2,swqrsm_1);
                swqrsm_3 = ltsamp(tt-3,swqrsm_2);
                swqrsm_4 = ltsamp(tt-4,swqrsm_3);
                if (swqrsm.lt_data - swqrsm_1.lt_data < onset && swqrsm_1.lt_data - swqrsm_2.lt_data < onset && swqrsm_2.lt_data - swqrsm_3.lt_data < onset && swqrsm_3.lt_data - swqrsm_4.lt_data < onset)
                    swqrsm = swqrsm_4;
                    tpq = tt - swqrsm.LP2n; %	/* account for phase shift */
                    break;
                end
            end

            if (learning~=1)
                % /* Check that we haven't reached the end of the record. */
                if tpq>length(data)
                    break;
                end
                % /* Record an annotation at the QRS onset */
                qrs = [qrs tpq];

                % J-points processing
                if (jflag)
                    tj = t+5;
                    for tt=t:t + fix(EyeClosing/2)
                        swqrsm = ltsamp(tt, swqrsm);
                        if swqrsm.lt_data > maxd - fix(maxd/10)
                            tj = tt;
                            break;
                        end
                    end
                    if tj>length(data)
                        break;
                    end
                    % Record an annotation at the J-point
                    jpoints = [jpoints tj];
                end
            end
            % /* Adjust thresholds */
            Ta = Ta + (maxd - Ta)/10;
            T1 = Ta / 3;

            % /* Lock out further detections during the eye-closing period */
            t = t + EyeClosing;
        end

    elseif (learning~=1)
    % 	    /* Once past the learning period, decrease threshold if no QRS
    % 	       was detected recently. */
        timer_d = timer_d + 1;
        if timer_d > ExpectPeriod && Ta > Tm
            Ta = Ta -1;
            T1 = Ta / 3;
        end
    end
    t = t+1;
end

end

function swqrsm = ltsamp(t,swqrsm)
    % /* ltsamp() returns a sample of the length transform of the input at time t.
    %    Since this program analyzes only one signal, ltsamp() does not have an
    %    input argument for specifying a signal number; rather, it always filters
    %    and returns samples from the signal designated by the global variable
    %    'sig'.  The caller must never "rewind" by more than swqrsm.BUFLN samples (the
    %    length of ltsamp()'s buffers). */
    while (t > swqrsm.lt_tt)
        swqrsm.Yn2 = swqrsm.Yn1;
        swqrsm.Yn1 = swqrsm.Yn;
        v0=swqrsm.data(1);
        v1=swqrsm.data(1);
        v2=swqrsm.data(1);
        if swqrsm.lt_tt>0 && swqrsm.lt_tt<=length(swqrsm.data)
            v0 = swqrsm.data(swqrsm.lt_tt);
        end
        if swqrsm.lt_tt-swqrsm.LPn>0 && (swqrsm.lt_tt-swqrsm.LPn)<=length(swqrsm.data)
            v1 = swqrsm.data(swqrsm.lt_tt-swqrsm.LPn);
        end
        if swqrsm.lt_tt-swqrsm.LP2n>0 && (swqrsm.lt_tt-swqrsm.LP2n)<=length(swqrsm.data)
            v2 = swqrsm.data(swqrsm.lt_tt-swqrsm.LP2n);
        end
        if v0~=-32768 && v1~=-32768 && v2~=-32768
            swqrsm.Yn = 2*swqrsm.Yn1 - swqrsm.Yn2 + v0 - 2*v1 + v2;
        end
        dy = fix((swqrsm.Yn - swqrsm.Yn1) / swqrsm.LP2n);	%	/* lowpass derivative of input */
        swqrsm.lt_tt=swqrsm.lt_tt+1;
        swqrsm.et = fix(sqrt(swqrsm.lfsc +dy*dy)); % /* length transform */
        id = mod(swqrsm.lt_tt,swqrsm.BUFLN);
        if id == 0
            id = swqrsm.BUFLN;
        end
        swqrsm.ebuf(id) = swqrsm.et;
        id2 = mod(swqrsm.lt_tt-swqrsm.LTwindow,swqrsm.BUFLN);
        if id2 == 0
            id2 = swqrsm.BUFLN;
        end
        swqrsm.aet = swqrsm.aet + (swqrsm.et - swqrsm.ebuf(id2));
        swqrsm.lbuf(id) = swqrsm.aet;
    % 	/* swqrsm.lbuf contains the average of the length-transformed samples over
    % 	   the interval from tt-swqrsm.LTwindow+1 to tt */
    end

    id3 = mod(t,swqrsm.BUFLN);
    if id3 == 0
        id3 = swqrsm.BUFLN;
    end
    swqrsm.lt_data = swqrsm.lbuf(id3);
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function g = glmderiv(net, x)
%GLMDERIV Evaluate derivatives of GLM outputs with respect to weights.
%
%	Description
%	G = GLMDERIV(NET, X) takes a network data structure NET and a matrix
%	of input vectors X and returns a three-index matrix mat{g} whose  I,
%	J, K element contains the derivative of network output K with respect
%	to weight or bias parameter J for input pattern I. The ordering of
%	the weight and bias parameters is defined by GLMUNPAK.
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Check arguments for consistency
errstring = consist(net, 'glm', x);
if ~isempty(errstring)
    error(errstring);
end

ndata = size(x, 1);
if isfield(net, 'mask')
  nwts = size(find(net.mask), 1);
  mask_array = logical(net.mask)*ones(1, net.nout);
else
  nwts = net.nwts;
end
g = zeros(ndata, nwts, net.nout);

temp = zeros(net.nwts, net.nout);
for n = 1:ndata
    % Weight matrix w1
    temp(1:(net.nin*net.nout), :) = kron(eye(net.nout), (x(n, :))');
    % Bias term b1
    temp(net.nin*net.nout+1:end, :) = eye(net.nout);
    if isfield(net, 'mask')
	g(n, :, :) = reshape(temp(find(mask_array)), nwts, net.nout);
    else
	g(n, :, :) = temp;
    end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [V,C] = calcVoronoi(v,varargin)
% compute the area of the Voronoi decomposition
%
% Input
%  v - @vector3d
%
% Output
%  V - list of Voronoi--Vertices
%  C - cell array of Voronoi--Vertices per generator
%
% See also
% voronoin

n = length(v);
v = reshape(v,[],1);

[x,y,z] = double(v);
faces = convhulln([x(:) y(:) z(:)],{'Qt','Pp','QJ'}); % delauny triangulation on sphere

% voronoi-vertices
V = normalize(cross(v.subSet(faces(:,3))-v.subSet(faces(:,1)),...
  v.subSet(faces(:,2))-v.subSet(faces(:,1))));

% voronoi-vertices around generators
[center, vertices] = sort(faces(:));

v = v.subSet(center);
vert = repmat(V,3,1);
vert = vert.subSet(vertices);

% the azimuth of a voronoi-vertex relativ to its generator
[ignore,azimuth] = polar(hr2quat(v,zvector).*cross(v,vert));

% sort the vertices clockwise around with respect to its center
[ignore,left] = sortrows([center azimuth]); %#ok<*ASGLU>

left = mod(vertices(left)'-1,length(V))+1;

% now we delete duplicated voronoi vertices
eps = 10^-10; % machine precision
[ignore,first,ind] = unique(round(squeeze(double(V))/eps)/eps,'rows');
V = V.subSet(first);
left = ind(left)';

% erase duplicated vertices in the pointer list
dublicated = find([diff(left)==0,false]);

% check whether the duplicated is in the next cell // they shouldn't be
% deleted
last = [0;find(diff(center));length(center)];
dublicated(ismember(dublicated,last)) = [];

left(dublicated) = [];
center(dublicated) = [];

C = cell(n,1);
last = [0;find(diff(center));length(center)];
for k=1:numel(last)-1
  ndx = last(k)+1:last(k+1);
  C{center(ndx(1))} = left( ndx );
end

author: mtex-toolbox

repo: mtex

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


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% Symmetry

% 	Even symmetry

syms x t k T
w=2*pi/T;
c1=(2/T)*int(x*sin(k*w*t),t,-T/2,0);

c2=(2/T)*int(x*sin(k*w*t),t,0,T/2);
c=c1+c2

t0=-2;
T=4;
w=2*pi/T;
syms t
x=t^2;
k=-5:5;
a=(1/T)*int(x*exp(-j*k*w*t),t,t0,t0+T);
stem(k,eval(a))
legend('a_k')

% 	Odd symmetry
figure
syms x t k T
w=2*pi/T;
b1=(2/T)*int(-x*cos(k*w*t),t,-T/2,0);

b2=(2/T)*int(x*cos(k*w*t),t,0,T/2);
b=b1+b2

t0=-2;
T=4;
w=2*pi/T;
syms t
x=t;
k=-5:5;
a=(1/T)*int(x*exp(-j*k*w*t),t,t0,t0+T);

stem(k,imag(eval(a)))
legend('a_k')


text:
%% UNIX_HEAD_TEST tests the UNIX facility for issuing operating system commands.
%
%  Discussion:
%
%    HEAD is a legal command on UNIX systems, which returns the first
%    lines of a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'UNIX_HEAD_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate a simple use of the "unix" command,\n' );
  fprintf ( 1, '  which allows MATLAB to invoke UNIX operating system commands.\n' );

  [ status, string ] = unix ( 'head unix_head_test.m' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Command status was %d\n', status );
  fprintf ( 1, '  Command output was %s\n', string );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Repeat the command, but add the -5 option.\n' );

  [ status, string ] = unix ( 'head -5 unix_head_test.m' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Command status was %d\n', status );
  fprintf ( 1, '  Command output was %s\n', string );

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


text:
function headingerr = getheadingerr(heading1,heading2)
%UNTITLED The input paras should in RAD
%   The input paras should in RAD
headingerr = RAD2DEG(heading1) - RAD2DEG(heading2);

headingerr(headingerr>180) = headingerr(headingerr>180)-360;
headingerr(headingerr<-180) = headingerr(headingerr<-180)+360;
end

repo: Wheel-SLAM


text:
function test225 ( )

%*****************************************************************************80
%
%% TEST225 tests L4MAT_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2011
%
%  Author:
%
%    John Burkardt
%
  m = 20;
  n = 50;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST225\n' );
  fprintf ( 1, '  L4MAT_PRINT prints a logical matrix.\n' );

  a = zeros ( m, n );

  for i = 1 : m
    for j = 1 : n
      a(i,j) = ( mod ( i, j ) == 0 );
    end
  end

  l4mat_print ( m, n, a, '  A(I,J) = I is divisible by J' );

  return
end


text:
function D = compute_pairwise_distance(X, Y)
    % efficient implementation taken from homework assignment 7
    if( ~isa(X,'double') || ~isa(Y,'double'))
        error( 'Inputs must be of type double'); end;
    m = size(X,1); n = size(Y,1);
    Yt = Y';
    XX = sum(X.*X,2);
    YY = sum(Yt.*Yt,1);
    D = XX(:,ones(1,n)) + YY(ones(1,m),:) - 2*X*Yt;
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function  commit()
%Commit Commit the currect directory to TortoiseSVN
% See also: svn
svn('commit',pwd)